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Abstract. 

Using a series of high-resolution N-body hydrodynamical numerical simulations, we investigate several scenarios 
for the evolution of the baryon budget in galactic halos. We derive individual halo star formation history (SFH), 
as well as the global star formation rate in the universe. We develop a simple analytical model that allows us 
to compute surprisingly accurate predicti ons, when compared to our simulations, but also to other simulations 
presented in ISpringel fc Hernauisd <2003bft . The model depends on two main parameters: the star formation time 
scale t* and the wind efficiency ??w We also compute, for halos of a given mass, the baryon fraction in each of 
the following phases: cold discs gas, hot halo gas and stars. Here again, our analytical model predictions are in 
good agreement with simulation results, if one correctly takes into account finite resolution effect. We compare 
predictions of our analytical model to several observational constraints, and conclude that a very narrow range 
of the model parameters is allowed. The important role played by galactic winds is outlined, as well as a possible 
'superwind' scenario in groups and clusters. The 'anti-hierarchical' behavior of observed SFH is well reproduced 
by our best model with t» = 3 Gyr and rj v = 1.5. We obtain in this case a present-day cosmic baryon budget of 
ft* ~ 0.004, ft coM ~ 0.0004, fihot ^ 0.01 and ft bac k ^ 0.02 (diffuse background). 
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1. Introduction 

Thanks to recent advances in the observation of high- 
redshift star-forming galaxies, and to the improved statis- 
tics in low-redshift galaxy surveys, it is now possible 
to have a quantitative view of the star formation his- 
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In the standard picture of galaxy formation, small 
mass objects collapse first, when gravitational instability 
enters the non-linear regime. This so-called 'hierarchical 
scenario' relies heavily on the observation of primordial 
density fluct uations imprinted on the Cosmic Microwave 
Background l|Spergel et al.ll20o"^) . These first generation 
objects eventually form stars, the so-called Population III 
stars, whose properties are likely to differ strongly from 
present-day galactic stars. The contribution of Population 
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III stars to the global star formation history in the uni- 
verse is currently a matter of debate: they likely con- 
tribute to a prompt initial enrichment of the intergalactic 
medium, as well as a possible early reionization epoch, 
but they may also cause a strong negative feedback on 
the star formation efficiency at this very early epoch. In 
this paper, we only adress star formation occurring inside 
galactic discs, concentrating on quiescent modes of star 
formation. This means that we only consider atomic pro- 
cesses to drive the cooling and subsequent catastrophic 
collapse of halo gas that ultimately leads to the forma- 
tion of a centrifugally supported disc. Molecular processes 
that are relevant for first stars formation, or for molecular 
clouds formation inside galactic discs are neglected in this 
study. 

Our main concern is therefore to follow the forma- 
tion of the first generation of gas rich, rotating discs, and 
the subsequent merging hierarchy that leads to larger and 
larger discs, up to present day spiral galaxies, as well as 
galaxy groups and clusters. The standard approach in cur- 
rent galaxy formation scenario is to consider that hot gas 
virializes first into extended dark matter halo potential 
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well. Depending on the Virial temperature of this hot gas, 
it cools rather rapidly, looses pressure support and col- 
lapses up to a point where centrifugal equilibrium sets in. 
It is inside these cold, rotating discs that star formation 
takes place at the end of a complex cascade of turbulent 
fragmentation and molecular cloud formation. On cosmo- 
logical scale, models of star formation still rely on heuristic 
recipes, partially based on observational constraints, and 
partially based on theoretical arguments. 

This rather simple picture is altered by feedback pro- 
cesses due to these stars. First of all, massive stars create 
collectively a strong UV background that photo-ionizes 
the intergalactic medium, preventing small mass halos 
from building their virialized, hot gas component, and al- 
tering their cooling efficiency. This means that, after the 
universe is re-ionized, we have to introduce a mass thresh- 
old (or equivalently a Virial temperature threshold) below 
which star formation is suppressed. This threshold allows 
us to distinguish between a diffuse component, defined by 
small mass dark matter halos with a very low gas fraction, 
and a star forming component, defined by high mass dark 
matter halos where disc and star formation is possible. 
The diffuse component is often called 'the Lyman alpha 
forest' or 'the smooth baryon background' in the litera- 
ture. In the present paper, we name the other component 
'the star forming halos'. Following standard definitions in 
galaxy formation studies, we further divide even more the 
baryons sitting in these 'galactic halos' into 3 different 
phases: hot gas, cold discs and stars. The purpose of this 
paper is to compute as precisely as possible the evolution 
of these 4 phases, for the universe as a whole, as well as 
for individual dark matter halos, thus studying the baryon 
logistics. 

A lot of other different feedback processes are currently 
under close examination in the literature, like supernovae- 
driven winds, young protostellar jets, AGN-driven jets 
and associated buoyant bubbles, and so on and so forth. 
Suc h outflows are actually observed in high redshift galax- 
ies They consist of high- velocity fountains 
of ionized gas, expelled from a parent, star-forming disc. 

These winds are likely t o be caused by collective ou tbreaks 

J a — • n II k 

of supernovae bubbles (de Avillez & Berry 2001). Other 

energy sources are however possible, like a central massive 
black hole, or collective jet ram pressure from star forming 
clouds. Regardless of their physical origin, these outflows 
are fundamental in explaining the metal enrichment in the 
intergalactic medium, as well as the high-metallicity ob- 
served in rich galaxy cluster. The modeling of such winds 
will also be addressed in this paper, in a simplified way, 
and their impact on the star formation history in the uni- 
verse as a whole will be computed. We will investigate in 
more detail the impact of winds on the star formation his- 
tory of individual star forming halos in a follow-up paper. 

Several approaches have been used to investigate the 
evolution of the baryons in the hierarch ical model of 
galax y f ormation. Numeric al simul ati ons llCen fe Ostrike: 
1992 : INavarro & White! Il993l iMihos fc Herriauis 
m& lKa,tzetal.l ll 99li lOnedin fc Ostrikerl Il997t 



Thacker fc Couchmanl [2000; 



Springel fc Hernauistl 



199 



199 



2003al) and semi-analy t ic mo d els llWhite fc Frenkl 
Somerville fc Primackl Il999t iKauffmann et alj 

:)le et al.ll2000t lHatton et alJl2003al) have been the most 
popular techniques. In a recent paper, an analytical 
appro ach has been proposed by iHernauist fc Springe! 
(2003) to compute the star formation history of the 
universe. Here, we present also a simple self-consistent 
analytical model. This model, validated on simulations, 
allows us to quickly compute the evolution of the 4 
baryons components, for the Universe as a whole, and 
also for an individual average halo. 

The paper is organized as follows: in the next section, 
we present numerical simulations of star formation in a 
ACDM universe, based on our Adaptive Mesh Refinement 
code named RAMSES. We then describe a simple analyti- 
cal model that allows us to compute the star formation his- 
tory inside individual halos and in the universe as a whole. 
In the section 4, we compare our model to simulation re- 
sults, showing that, once calibrated to our simulations, the 
model works very well in predicting the evolution of the 4 
different baryon phases. We finally compare our model to 
several observational constraints. Our analytical approach 
allows us to explore efficiently the parameter space, which 
is in our case limited to a 2-dimensional space: star forma- 
tion time i* versus wind efficiency r/ w . We finally confront 
our predictions to current observational constraints, and 
found a rather narrow parameter range. 



2. Simulations 

Our simulations were performed using the Adaptive 
Mesh Refin ement co d e cal led RAMSES and described 
in detail by iTevssierl (2002). The N-body solver is very 
similar to the ART code l|Kravtsov et al ] Il997h and 
the hydrodynamical solver is based on a second-order 
Goduno v-type meth od, called the MUSCL-Hancock 
scheme dTorol 119971) . We evolve the collisionless dark 
matter particle distribution by solving the Vlasov-Poisson 
equation, and the baryonic component by solving the 
Euler equation with gravity source terms. More detail 
about our hydrodynamical solver and our refinement 
strategy are given in appendix ^] We also solve for a 
heating and cooling source terms in the energy equation, 
ass uming primordial H and He plasma photo-ionized by 
the lHaardt fc MadaiJ l)l99fil) UV background. In dense 
and cold regions of the flow, we turn a fraction of the 
gas into collisionless 'star' particles. This numerical 
approach is widely used in current galaxy forma- 
tion st udies llOen fc Ostrikerl Il992t iNa'varro fc Whitf " 



tion st udies llOen fc Ostrikerl Il992t INavarro fc 
1995 IMihos fc Hernnuisd Il994l lKa,tz et alJ 
Gnedin fc Ostrikerl Il997t iThacker fc Couchmanl 



199 
200C 



We will briefly recall here the properties of our own im- 
plementation. We then describe cosmological parameters, 
box sizes and mass resolution we use in our three main 
simulation series. 
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Name 


L10N64S30 


L10N128S30 


L10N256S30 


L10N512S30 


L(/i- i Mpc) 


10 


10 


10 


10 




0.036 


0.036 


0.036 


0.036 


to(Gyr) 


30 


30 


30 


30 




3 


3 


3 


3 


A^ccll 


64^ 


128^ 


256 a 


512 a 


p 

<-max 


5 


5 


5 


5 


^end 


3 


3 


3 


3 


m M(h~ L M @ ) 


2.8 x 10 s 


3.4 x 10 Y 


4.3 x 10 b 


5.4 x 10 b 


Ax(h~ L kpc) 


4.9 


2.4 


1.2 


0.6 



Table 1. Runtime parameters for the 'convergence study' simulation suite. L is the box length, no is the density threshold at 
z = 0, to is the star formation time scale for density no, ao gives the evolution of the star formation time scale with redshift, 
iV ce ii is the number of cells at the coarse level, ^max is the level of refinement, z cn d is the final redshift, tudm is the dark matter 
particle mass and A x is the spatial resolution. 



2.1. Star formation recipe 

We now describe our method for implementing star for- 
mation in the RAMSES cosmological code. It is based on 
an heuristic approach adopted in many, if not all, cosmo- 
logical studies. For a complete review of various imple- 
mentations, see lKav et alJ l)2002l) . Basically, one considers 
that star formation proceeds at a given time scale, written 
here t* , in region where one or several physical criteria are 
fullfiled. In this paper, as well as in many previous other 
papers in the literature, we adopt a simple scheme to turn 
gas mass into star particles, adding a source term in the 
continuity equation 



Dp 
Dt 
Dp 
Dt 



= ifp>p (i + ^r\ 

= otherwise, 



(1) 



where the threshold density may depend on redshift 
through index ao. The star formation time scale t* is pro- 
portional to the local free-fall time 



t* — to ( — 
Pa 



-1/2 



(2) 



In the literature, one can find basically two different ap- 
proaches: a density threshold constant in physical units, 
corresponding here to ao = 0, and a density threshold con- 
stant in comoving units, corresponding here to ao = 3. 

This simple model of star formation has bee n discussed 
and criticized extensively in the literature IjKav et alJ 
l2002h : we recall here briefly its possible physical and ob- 
servational origins. Spiral galaxies in our nearby universe 
are seen to form sta rs at a rate given by the Kennicutt 
law l|Kennicuttlll998j) . similar to equation ^ with volume 
densities p and po replaced by average disc surface densi- 
ties. 



S = 10 M /pc 2 and t = 2 Gyr. 



(3) 



The physical origin of such a behavior is not clearly 
identified yet. One good candidate is the sustained, in- 



terstellar, self-gravitating turbulent cascade, which con- 
trols the mass flux between large scale filaments in the 
disc and small scale star forming molecular clouds, at a 
rate given by the loca l (on large scale, though) free-fall 
time (|Elmegreenll2f)f)2|l . Such a precise description of star 
formation is completely irrelevant for our present cosmo- 
logical study. Using a 'sub-cell approach', similar in spirit 
to the one developed by the fluid mechanics community 
and often named the k — e method, our current modeling 
of star formation try to mimic in a statistical sense the 
complex behavior of the in terstellar medium. 

Alon g the same lines. lYepes et all l)l997ft and more 
recently ISpringel fc Hernauistl l|2003aj) have developed a 
multiphase model of the interstellar medium, based on 
iMcKee fc Ostrikerl 1 1977f) early work on molecular clouds 
evaporation within supernovae remnants hot bubbles. 
This multiphase model offers an interesting alternative to 
the previous turbulent model. It also gives the possibility 
to computes self-consistently the star formation parame- 
ters, which are 



rip = 0.1 cm and ao = and to — 2.1 Gyr. 



(4) 



Recall that in this case, star formation is allowed in 
region whose gas density lies above a physical density 
threshold. In cosmological simulations, however, the co- 
moving density is usually preferred to define collapsed, 
high-de nsity clumps where star form ation is likely to oc- 
cur. In ISpringel fc Hernauistl l|200 3b) for example, equa- 
tion ^ is augmented with the requirement that the local 
overdensity should exceed 200. This guarantees that star 
formation cannot occur in smooth regions of the cosmo- 
logical flow, but only within collapsed, virialized halos. 

In the present paper, we would like to explore a dif- 
ferent o ption than the one used in Springcl & Hcrnauist 
( 20033). We use the following star formation density 
threshold 



no = 0.036 cm and ao = 3. 



(5) 



This corresponds to a baryon overdensity threshold of 
1.6 x 10 5 . This approach was also adopted in earlier 
works on the star formation history in the universe 
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llCen fc Ostrikerl Il992t iKav et all [20021: iNagamine et alJ 
120041) . with threshold overdensities ranging from 5 to 10 , 
depending on the authors. In our case, this is motivated 
by our AMR refinement scheme: n Q (l + z) 3 corresponds 
exactly to the density threshold triggering the maximum 
level of refinement £ max = 5. When baryons cool and set- 
tle down at the centre of their host dark matter halos, 
the density does increase dramatically. Our AMR code 
describes this collapse accurately, by adding recursively 
new cells at the centre of the halo. At the point when the 
maximum level of refinement has been reached, the numer- 
ical description of the collapse is not valid anymore. We 
then turn on star formation, as our sub-cell modelling of 
this uncomplete collapse, that would have lead ultimately 
to the formation of one or several star forming molecular 
clouds. 

As soon as star formation is active, we create collision- 
less star particles of constant mass m*. As we explain in 
appendix IA. 21 this constant mass is chosen to be equal to 
the initial mass resolution in the gas distribution, in order 
to prevent spurious refinement or de-refinement triggered 
by star formation. Using a constant mass for our star par- 
ticles has also the advantage of controlling the maximum 
number of star particles created at the end of the sim- 
ulation. In order to solve for the star formation source 
term (Eq.QJ with constant star particles mass, we need to 
a dopt a stoch astic approach, similar to the one proposed 
in lKatd l|l992F) . In star forming cells, we generate N equal 
mass star particles, where N is drawn from a Poisson pro- 
cess, with probability 



P(7V) = — exp(-A), 

and with parameter (or mean value) 
/ pAx 3 \ At 



A = 



V m * 



(6) 



(7) 



These star particles are created at each time step, at the 
very centre of their parent cell. They are given a velocity 
equal to the local fluid velocity, plus a random component 
that we take equal to the local sound speed in the gas. 
The corresponding mass, momentum and internal energy 
is of course removed from the parent cell's conservatively. 

In the simulation presented here, the only free param- 
eter is the star formation time scale at comoving threshold 
density, to- 

2.2. Parameters 

We assume throughout this paper a single cosmological 
model, the so called 'concordance model', with a cold dark 
matter component with Q m = 0.3, a baryon component 
with fib = 0.04 and a dark energy component with 17a = 
0.7. The Hubble constant was set to h = 0.7. The initial 
power spectrum is assumed to be Harrisson-Zeldovich with 
n = 1 and a normalization constrained by as — 0.93. The 



exact fu nctional form we use for the transfer function is 
given in ISueivamal l|l995() . 

We performed several AMR simulations, varying the 
box size, the number of particles and the star formation 
time scale to- The density threshold for star formation was 
set to no = 0.036 cm -3 , and was held constant in comoving 
units. In our notation, this translates into ctQ = 3. This 
star formation threshold corresponds to the overdensity 
threshold triggering our last level of refinement ^ max = 5. 

Our simulation suite can be organized in 3 different 
sets. The first set was designed to perform a convergence 
study. We used a box size of 10/i -1 Mpc, with identical 
physical and cosmological parameters. We vary only the 
initial number of particles and grid cells, from 64 3 to 512 3 . 
Runtime parameters for this 'convergence study' simula- 
tions are summarized in Tabled 

The second set of simulations was designed to explore 
a 'high efficiency' SFR scenario with the maximum dy- 
namical range we can afford (Table EJ. We used for that 
purpose 3 different box size: 1, 10 and 100 h _1 Mpc. All 
these simulations were performed using 256 3 particles and 
the same number of grid cells. This translates into a dark 
matter particle mass of 4.3 x 10 3 , 4.3 x 10 6 and 4.3 x 10 9 
h _1 M0 . The third set of simulations is similar to the sec- 
ond one but for a 'low efficiency' SFR scenario (Table 0J). 

The 2 last simulation sets span a huge dynamical range 
in mass, but each simulation is valid for a limited range of 
redshifts. When the typical scale of non linearity reaches 
the box size, the simulation has to be stopped: the sim- 
ulated volume is not statistically representative of the 
universe as a whole, and spurious effects due to periodic 
boundary conditions become visible. The stopping redshift 
was chosen to be 5, 2.5 and 0, for a box size of 1, 10 and 
100 h _1 Mpc respectively. 

We have not implemented feedback processes in the 
RAMSES code: this is currently under development. As 
explained in the introduction, galactic winds are a key 
ingredient in computing the star formation history in the 
universe. In this paper, we use SPH sim ulations results 
obtained bv lSpringel fc Hernauistl l|2003bh to estimate the 
influence of winds in the overall baryon budget. 

Three set of AMR simulations performed by us us- 
ing the RAMSES code ('convergence study', 'high effi- 
ciency ' and 'low efficiency') and o ne set of SPH simula- 
tions l)Springel fc Hernauistll200 3b'l compose the data we 
analyze and discuss intensively in this paper. 



2.3. Results 

We present now general results obtained by the RAMSES 
code, for a typical cooling and star formation run. We 
would like to outline that current high-resolution numeri- 
cal simulations reproduce qualitatively the global picture 
of galaxy formation: fast cooling gas builds up centrifu- 
gally supported discs at the center of dark matter haloes, 
in which star formation quietly proceeds. We will analyze 
our results more quantitatively in the next sections. 
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Name 


L100N256S3 


L10N256S3 


L1N256S3 


L(/i- i Mpc) 


100 


10 


1 




0.036 


0.036 


0.036 


to(Gyr) 


3 


3 


3 




3 


3 


3 


A^ccll 


256 a 


256 a 


256 a 


p 

<-max 


5 


5 


5 


^end 





2.5 


5.5 


m D M(/i _1 M ) 


4.3 x 10 a 


4.3 x 10 b 


4.3 x 10 d 


Ax(h~ L kpc) 


12 


1.2 


0.12 



Table 2. Runtime parameters for the 'high efficiency' simulation suite. The meaning of the symbols is the same as tabled 



Name 


L100N256S30 


L10N256S30 


L1N256S30 


L(ft _1 Mpc) 


100 


10 


1 


no(cm _il ) 


0.036 


0.036 


0.036 


to(Gyr) 


30 


30 


30 


QO 


3 


3 


3 


iVccll 


256 a 


256 a 


256 a 


^max 


5 


5 


5 


Zend 





2.5 


5.5 


moM(h~ L M e ) 


4.3 x 10 y 


4.3 x 10 b 


4.3 x 10 a 


Ax(/!, _i kpc) 


12 


1.2 


0.12 



Table 3. Runtime parameters for the 'low efficiency' simulation suite. The meaning of the symbols is the same as tabled 
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Fig. 1. Gas density (up) and stellar density (down) in our highest resolution run L10N512S30. The first map on the left of each 
row shows a projection of the full periodic volume. The squares in each map delimits the zoom region where the next image in 
the row is defined. For color map definitions, see text. 



Figure n presents the simulated density field projected 
along one principal axis of the comoving periodic box. Run 
L10N512S30 is shown, at redshift z = 2.8. The projected 
gas density is shown with a logarithmically scaled color 
table. The 4 images show a zooming sequence, starting 
from the whole periodic box with the typical large scale 



filamentary structure, down to a particular region where 
2 spiral discs are clearly visible. 

The second set of images in Figure ^ shows the same 
zooming sequence for the projected stellar density. These 
images were computed using a 'true color' color scheme: 
stars are divided into 3 families according to the formation 
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redshift of each star particle (3 < Zf < 4, 4 < Zf < 5 
and 5 < Z{). 3 projected density maps are computed for 
each family. The 3 images are then combined using the 
RGB color scheme: 'Red' stands here for old stars, 'Blue' 
for young stars and 'Green' for intermediate formation 
redshift stars. 

These rather spectacular images compare favorably 
with observed high redshift galaxies observed for example 
with the Hubble Space Telescope. We have to be careful in 
making definitive conclusions: the spiral discs we observe 
in our simulations arc highly underresolved. Figurenillus- 
trates this point: the AMR grid structure clearly shows up 
in the highest resolution figure, demonstrating that these 
runs are not meant to resolve the internal structure of 
galactic discs. We are however confident in the computed 
baryon fraction inside each halo, as soon as the halo mass 
is greater than a few hundred particles. 

The traditional way of describing the baryon density 
field is to divide it into 4 phases. This phase separa- 
tion is apparent in Figure El which shows the density- 
temperature histogram at z = for our large box run 
(100 h" 1 Mpc) L100N256S30. The first regime occurs at 
low density (p < 200p) and low temperature (T < 10 5 
K) in the phase space diagram. This is the diffuse inter- 
galactic medium, also known as the Lyman alpha forest. 
The diffuse intergalactic UV flux is responsible for pre- 
venting this warm gas from collapsing into their parent 
dark matter halo. 

The second regime occurs at high density (p > 10 5 p) 
and low temperature (T < 10 5 K), and corresponds in 
our simulations to cold, centrifugally supported discs. This 
gas is termed here 'cold gas'. It is shown in Figure as 
the rightermost box. In this region, the quasi-neutral gas 
follows a very tight p-T relation, typical of high-density 
HI discs. 

The remaining gas corresponds to shock heated gas 
into virialized halos. We call it 'hot gas' although, when 
cooling is efficient, its temperature never exceed a few 10 
K. This gas, in rough hydrostatic equilibrium, spans a 
large density range, from p ~ p in the outskirts of dark 
matter halos up to p ~ 10 5 p in the X ray emitting cores. 
This ionized gas is also rapidly cooling from the inside out, 
leading to the formation of the cold phase. 

The last component in the baryon budget is the stel- 
lar phase. Stars originate from the highest density tail in 
the phase space diagram. In our model, star formation 
occurs above the overdensity threshold p > 10 5 p, which 
corresponds to our cold phase. 

Molecular cooling as well as supernovae feedback are 
not modeled here. The cold neutral gas is in reality de- 
composed into several components, typical of the inter- 
stellar medium: molecular, very cold clouds embedded 
into hot, supernovae driven bubbles. The correct descrip- 
tion of this multiphase interstellar medium is far beyond 
the possibilities of current cosmological simulations. The 
only possibility is to rely on sub-c ell modeling , alonp 
the lines described for ep 
and lSpringel fc Hernauistl i 



. These authors noticed 



that their multiphase modeling does not affect the over- 
all baryon budget between the hot halo gas, the cold disc 
component and the stellar fraction. The main effect of this 
multiphase approach was to modify the internal structure 
of the gaseous discs, altering the effective equation of state 
of the dense gaseous component. 

The most important feature that is likely to affect the 
computed baryon budget is the mass resolution of the 
code. We have to carefully assess its effect on our results. 
Let us first examine Figure |3 which shows a sequence 
of images of the projected gas density fields in a chosen 
region of the 'convergence study' suite. Initial conditions 
were generated self-consistently for these 4 simulations, in 
order to recover the same large scale distribution. From 
the lowest resolution run with 64 3 particles to start with, 
up to the highest resolution run with 512 3 dark matter 
particles, one clearly sees the spectacular increase in small 
mass haloes along the filaments and in the voids in be- 
tween. 

If we now examine in Figure the corresponding stel- 
lar density distributions, we see that small haloes appear 
devoid of stars. Small haloes are indeed unsufficently re- 
solved to reach the high density contrast required to allow 
star formation. Cooling is also very likely affected by the 
poor resolution in these small halos. We will see in the 
next sections that for our simulations the minimum mass 
for a halo to host stellar particles lies around 400 dark 
matter particles. 

Another effect of mass resolution is also visible in the 
same figure. Large galactic discs obtained at a given res- 
olution tend to fragment at higher resolution, leading to 
smaller discs with several orbiting satellites. These satel- 
lites are the remnants of progenitor halos, which form ear- 
lier at smaller mass. Insufficient resolution also affects the 
history of mass assembly inside galactic halos. Figure |31 
shows that the color of a given galaxy is affected by the 
finite mass resolution. At low resolution, star formation is 
artificially delayed to late times and the galaxy appears 
blue, while at higher resolution, the correct star formation 
history is recovered and the same galaxy appears red. The 
purpose of this paper is to carefully estimate the effect of 
finite mass resolution on our predictions. 

3. Model 

In this section, we present a simple analytical model to 
compute the evolution of the baryon budget in the uni- 
verse. The purpose of this analytical model is to shed 
light on the complex behavior of our numerical simula- 
tions. Our analytical treatment of the cosmic star forma- 
tion history is, of course, unaffected by finite resolution 
effect. It is however a crude model, and a careful compari- 
son with numerical simulations is required to validate our 
approach. 

This model differs with the 'semi-analytic' model- 
ing of galaxy formation, an approach pursued by several 
teams llWhite fc Frenklll99lt ISomerville & Primacklll999l 
iKauffmann et alJ Il999t ICole et alJ l2000t lHatton et all 
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Fig. 3. Projected density maps (top) and projected stellar density (bottom) of our 'convergence study' simulations series 
(L = 10 h _1 Mpc). The image size is 1.25 h -1 Mpc wide on a side. From left to right, the number of dark matter particles is 
increased from 64 3 , to 128 3 , 256 3 and 512 3 . 




I 



10" 5 



2 4 

log(p/<p» 

Fig. 2. The different baryon phases in the p — T diagram. 
Gray contours show a mass-weighted histogram: the baryon 
mass fraction at a given density and temperature. Each region 
corresponds to a given phase (diffuse background, hot or cold 
gas), as defined in the text. 



In this paper, our goal is to compute the evolution 
in the mass fraction of the 4 different components in 
the baryon distribution: diffuse background, hot virializcd 
plasma, cold neutral discs and stars. We will make predic- 
tions for the universe as a whole, but also for individual ha- 
los of a given mass. Semi-analytical models, when coupled 
to N body simulations, can predict the baryon history for 
individual halo, based on the specific merging history at 
the origin of the halo hierarchical mass assembly. Our sim- 
ple model allows us to compute only the average baryon 
history for halos of a given mass range. 

In a recent paper, an analytical approach has been pro- 
posed to compute the star forma tion history in the uni- 
verse l|Hernauist k Springe]| l2003) . The proposed method 
was to use the Press-Schechter theory for the dark matter 
halos statistics together with the star formation rate as 
a function of halo mass, as measured in SPH numerical 
simulations. In the present paper, we develop a fully self- 
consistent analytical model, slightly more c o mplex than 
the one proposed by iHernauist fc Springell l)2003|) . but 
based on a similar approach. This self-consistency allows 
us to compute the mass fraction evolution of the 4 baryon 
components, for the universe as a whole, and also for an 
individual average halo. 



2003a). These models are based on a quite sophisticated 
treatment of the physics of galaxy formation: cooling, 
star formation and spiral discs evolution are few exam- 
ples among the numerous ingredients in semi-analytical 
modeling. 



3.1. Halo Model 

The method we use in this paper to compute the star for- 
mation his_toryJs_l2a^ed OTi_what is usually called the halo 
model jCoorav fc Shetbll2002t . The idea is to decompose 
dark matter and baryon density fields into a collection 
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of virialized halos. whose distribution is described by the 
IPress fe Schechte mass function. This approach 

was introduced to describe galaxy and dark matter clus- 
tering, using two additional ingredients: linear halo bias- 
ing the ory and the Navarro, Frenk fc White (1995) density 
profile l|Ma fc Fr\l2000HSeliakl20QO|) . Later on, several au- 
thors used similar tools to compute the Sunyaev-Zeldovich 
power sp ectrum, as well as the mean Comptoniza tion pa- 
rameter lJCooravll200rl iRefregier fc Tevssierll2002l) . 

A feature common to most of these earlier works is 
the rather static picture they give to dark matter halos 
evolution. We propose here to use a similar approach, in 
order to compute the history of star formation within each 
dark matter halo. Our methodology differs somewhat from 
earlier works, since it is based on a two step approach. 

In the first step, using static halo model, we compute 
the mass transfer rates between each phase. In the sec- 
ond step, using these computed mass transfer rates, we 
solve for the mass fraction evolution equations, a system 
of first order Ordinary Differential Equations (ODE). Our 
analytical model is therefore based on the computation of 
the baryon logistics. In order to compute the baryon mass 
fraction locked up in stars for example, we need to solve 
the complete set of ODE from a large redshift (z = 200 
say) down to z=0. 

In this paper, we define the halo mass M as A^oo, 
the total mass enclosed in radius i?200i where the mean 
overdensity (relative to the background density) is A = 
200. 

M = M 200 - —p(z)AR 3 200 . (8) 



This choice differs from earlier definitions, where A was 
either defined relative to the critical density, or A was 
a function of redshift, as suggested by the sp heri cal col- 
lapse m odel. As shown bv lJenkins et alJ l)200l[) and I White! 
(2002), these earlier definitions are not suited to the Press 
& Schechter approach. Although they are based on physi- 
cal principles, they destroy the self-similarity of the Press 
& Schechter mass function. 

In this paper, we therefore adopt A = 200 (relative 
to the mean background density) independent of redshift. 
We also use very often the halo circular velocity V200 and 
the halo Virial temperature T200, instead of the halo mass. 
The Virial temperature must not be considered as a true 
physical temperature, but rather as yet another mass pa- 
rameterization. In this paper, the Virial temperature is 
defined as 

/xm H GM200 , , 

l 9 J 



R2 



(with fj, the mean molecular weight) and the circular ve- 
locity as 



V200 



GM 2 , 



R2 



(10) 



In the 5 following sections, we are going to compute 
the cosmic rates between the 4 phases. If necessary the 
reader can go directly to section 1X71 




10" 
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Fig. 4. Time evolution of the Minimal Mass M m in for two 
reionization scenarios: 2 r =20 (solid line) and z r =6 (dashed 
line). In both case, the background temperature after reion- 
ization was set to T r = 6 x 10 K. The Minimal Cooling tem- 
perature was also set to T coo i = 6 x 10 3 K. 

3.2. Minimal Mass 

The first component in the baryon budget is the diffuse 
background. This may be the most important one, since 
it is the reservoir of fresh gas that will eventually feed star 
forming halos at all epochs. It is usually called the Inter 
Galactic Medium (IGM) or the Lyman Alpha Forest. We 
need to give a precise definition of what we call 'diffuse 
background' in this paper. 

The diffuse background is the baryon component as- 
sociated to dark matter halo with masses lower than the 
minimal mass M m i n , below which cooling is inefficient and 
pressure forces prevent baryons to collapse in their poten- 
tial wells. This minimal mass is therefore fundamental be- 
cause it is the transition between 'star forming halos' and 
'diffuse' one. 

This Minimal Mass is taken to be the maximum be- 
tween the Filtering Mass and the Minimal Cooling Mass. 
The filtering mass Mp is the minimal halo mass above 
which baryons are able to fall into their dark matter halo 
potential wells (see Appendix [BJ) . The minimal cooling 
mass Af coo i is the mass above which the gas is able to cool 
and therefore to form stars (see Appendix Q. 

We implement this using a smooth function of both 
Virial temperatures 



T, 



cool- 



(11) 



The Minimal Mass is finally computed using the Virial 
relation (Eq. |5J). We plot in Figure 0] M m j n as a func- 
tion of redshift, for our two extreme reionization sce- 
narios. At early times, this mass remains roughly con- 
stant, independent of redshift. As reionization proceeds, 
this Minimal Mass increases steadily, up to a rather large 
value M min ~ 10 11 h _1 M Q today. 
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3.3. Cosmic Accretion Rate 

Using the Press & Schechter formalism, we now compute 
the mass fraction in the diffuse background and the mass 
fraction in star forming halos. Since the Minimal Mass 
is considered here as the mass threshold between these 
two components, and assuming that the baryon fraction 
in each halo is equal to the universal one, we have 



/hot = f(M> M min ) = / b jf ^ exp 



with 



S c (t) 
cr(M m in) 



and S c (t) 



1.686 

d+WY 



(-v 2 /2)dv, (12) 



(13) 



where D + is the linear growth factor and cr(M mill ) the 
variance of the density field smoothed at the Minimal 
Mass scale. The rate at which baryons are transferred from 
the diffuse background to star forming halos is computed 
by taking the time derivative of the previous equation 



/acc 



4fi 



hot 



4fi 



back 



dt 



dt 



= -/tA 



-exp(- I / in /2),(14) 



We then define the Cosmic Accretion Rate of fresh diffuse 
gas into star forming halos by 

/acc = W acc / back , (15) 

where the accretion rate, in units of Gyr , is given by 



2 exp(-iA n /2) 
7T erfc(^ min /\/2) 



(16) 



This last equations give the mass accretion rate of diffuse 
gas into star forming halos in the general case, for which 
the mass fraction in the diffuse component is allowed to 
vary from its canonical value. Note that this accretion rate 
has nothing in com mon to the traditiona l mass accretion 
rate on a given halo IjLacev fc C ole 1993). This rate gives 
the fraction of fresh diffuse gas dispatched among all star 
forming halos. This fresh gas in assumed to be transferred 
exclusively to the hot plasma component. The two vari- 
ables /back and /hot refer therefore to the total mass frac- 
tion in the background and the total mass fraction in the 
hot component, both integrated over the PS distribution. 

It is also possible to compute the Cosmic Accretion 
Rate on a halo b y halo basis, using the Extended Pres s 
Schechter theory (|Bond et al.lll99ltlLacev fc Coldll993|) . 
This theory allows to compute the progenitors mass dis- 
tribution as a function of time, for a given parent halo 
mass Mo, up to the 'halo formation time' to. The individ- 
ual Cosmic Accretion Rates are very similar to the one 
computed for the whole universe. We follow the same pro- 
cedure, computing first the mass fraction in star forming 
halos, assuming that each progenitor hosts a baryon frac- 
tion equal to the universal one. 



/hot (Ml, t ) — / b 



— exp(-j/ /2)dv, 



(17) 



10.000 




2.0 



Fig. 5. Cosmic Accretion Rate for the ACDM cosmology with 
z r = 20 for the universe as a whole (thick solid line) and for 
various halo masses (thin lines). Halo mass are, from top to 
bottom, Mo = 10 13 , 10 12 , 10 11 , 5 x 10 10 , 2.5 x 10 10 and 10 10 
h~ Mq. The halo formation redshift is set to zo = 0. 



with this time 



^min (M , t ) 



5 c (t)-5 c (t ) 



y/(r(M n 



cr(Mo) 



(18) 



The accretion rate is then computed exactly as for the 
previous case, using Equations 1151 and 1161 with however 
a different value for f m in given by Equation 1181 We want 
to stress again that we do not consider accretion of satel- 
lite halos on the most massive progenitor, which is the 
traditional way of computing the accretion rate. Here, we 
consider accretion of diffuse gas on all star forming pro- 
genitors of the final halo. 

This fresh gas contributes to fill up dark matter ha- 
los with hot, virialized gas. Hot gas coming from satellite 
halos is automatically accounted for in our formalism. If 
we assume very short cooling rates and instantaneous star 
formation, this Cosmic Accretion Rate is nothing but the 
Star Formation History in the universe. In a realistic case, 
star formation and cooling introduce a delay in the curve. 

This cosmic accretion rate depends on the thermal his- 
tory of the background gas, on the density power spectrum 
through u(M), and on the cosmological model through 
D + (t). Figure [5] shows the accretion rates (in Gyr -1 ) as a 
function of redshift, for a ACDM universe, and for differ- 
ent halo masses. The halo formation redshift was fixed to 
z Q = 0. 

For small mass halos, accretion stops abruptly as the 
Minimal Mass reaches the parent halo mass. This means 
that star forming progenitors are not present anymore in 
the halo. For halos more massive than M m j n , accretion 
remains active up to the halo formation redshift. The ac- 
cretion rate actually diverges as z — > zo, as all the remain- 
ing diffuse gas is accreted into the final virialized halo. 



10 



Yann Rasera and Romain Teyssier: The History of the Baryon Budget 



This diffuse mass accretion rate remains however finite 
for Mq > M m j n , and we can compute its value as z — > zq 



/a' 



-h 



Sc(t ) 



* ^(M min ) 2 - <r{M y 



(19) 



3.4. Cosmic Cooling Rate 



Our third baryon component (namely cold atomic gas in 
centrifugally supported discs) is progressively built up by 
accreting cooling gas into the very center of their parent 
dark matter halo. We need to estimate the global rate 
at which hot gas is transferred into this dense and cold 
component. Using EPS theory, we compute this rate for 
individual halo mass (Mo, Zq). Note that we recover the 
results for the Universe as a whole by taking the limit 
Mq — ► +oo and zq — » — 1. 

We follow our basic methodology, assuming this time 
that all baryons are in the hot halo phase. Using our sim- 
ple cooling model detailed in Appendix [U] we compute 
the total amount of gas cooling from the hot halo compo- 
nent during a unit time interval by integrating the instan- 
taneous cooling rate over the PS mass distribution from 
Af m i n to M max , the Maximal Cooling Mass. Indeed, above 
this mass (see Appendix 0, the halo enter in the slow 
regime cooling. We therefore neglect this contribution. It 
gives 



/cool — fb . 

^orb 



- exp(—u /2)dv, 



(20) 



where v m - m is defined by Equation 1181 f max corresponds 
to the Maximum Cooling Mass M max 



(M , t ) = 



5 c (t)-6 c (t ) 



(21) 



vMMnax) 2 " CT(Mo) 2 ' 

and t rb is the orbital decay timescale (see Appendix 0. 
We then define the Cosmic Cooling Rate of hot halo gas 
into cold gaseous discs by 



We therefore consider star formation in a dark matter 
halo as a function of the total amount of cold gas in that 
halo. The star formation rate in each halo is simply given 
by M„ = o^Mcoid with is the average star formation 
rate in that halo. In order to compute this average time 
scale from first principles, one needs to integrate the local, 
density dependent, star formation rate over the cold gas 
density PDF. 

In this analytical model, however, we consider star for- 
mation models inspired by star formation recipes used 
in numerical simulatio ns and by semi-analytical models 
i Somerville et al.ll200l[) . The halo star formation rate is 
parameterized by 



w* = i (l + zf^ 2 , 



(24) 



where t* is the present day star formation time scale and 
a* is the acceleration parameter. In the literature, two ba- 
sic quiescent models are usually discussed in galaxy for- 
mation studies. 

The first model, usually referred to as a 'constant effi- 
ciency' model, assumes that the halo star formation time 
is a constant a* = 0. This model corresponds to numerical 
simulations with a constant star formation density thresh- 
old. 

The second model assumes that a* = 3. It is usually 
called an 'accelerated efficiency' model. The star forma- 
tion time scale decreases with redshift (as the mean den- 
sity of the Universe increases). This model corresponds 
to numerical simulations with a constant star forma- 
tion overdensity threshold. It is also used in semi- ana- 
lytical models to mim ic starbursts triggered by mergers 
llSomerville et alJl200ll) . 

We compute now the global star formation rate, using 
our basic methodology. Since the halo star formation rate, 
in our simple scenario, does not depend on halo mass, we 
can integrate over the EPS mass function, and obtain the 
Cosmic Star Formation Rate as 



/cool — ^cool/hot, (22) /* — W*/cold- 

where the cooling rate, in units of Gyr , is given by 
I erfc^min/ y/2) - erfc ( v max / \f% 



(25) 



^cool 



erfc(f min /v / 2) 



(23) 



3.6. Cosmic Winds 



The Cosmic Cooling Rate depends on the cosmological 
model, on the thermal history of the background and on 
the details o f the cooling model. A similar model has been 
proposed bv lvan den Boschl l)2002f) . in a different context. 

3.5. Star formation models 

In this simple analytical model, we completely discard the 
description of the gaseous discs. Predicting the disc sizes 
and surface density profiles obtained in the hierarchical 
scenario of structure formation is beyond the scope of this 
paper. We are only interested in the global baryon budget, 
and more precisely in the global star formation history. 



The last ingredient in our model, but not the least, is 
the contribution of galactic winds to the overall baryon 
budget. It is a well known issue in current models of galaxy 
formation that without feedback processes, most baryons 
would end up into cold gas or stars, in contradiction with 
several observational constra ints. This problem is k nown 
as the 'overcooling problem' l|Blanchard et al.lll992f) . 

As discussed in the introduction, the exact nature of 
the dominant feedback process is still unknown. It is most 
likely that various processes are in competition, and their 
impact o n baryons may vary as a func tion of halo mass. 
Following ISnrineel fc Hernauistl l)2003bft . we assume m our 
model that winds occur during star formation events, 
probably related to supernovae. We therefore assume that 
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Fig. 6. Cosmic Accretion Rate (solid line), Cosmic Cooling 
Rate (dotted line), Cosmic Star Formation Rate (dashed line) 
and Cosmic Outflow Rate (unbound fraction, dot-dashed line) 
for the ACDM cosmology with the following model parameters: 
z r = 20, U = 3 Gyr, a* = 0, T w = 2 x 10 6 K and r) v = 1.5. 
These rates were computed using Afo = +oo and zq — — 1, and 
therefore corresponds to the universe as a whole. 



Fig. 7. History of the mass fraction in the diffuse background 
(solid line), in the hot halo gas (dotted line), in the cold discs 
(dashed line) and in stars (dot-dashed line) for the ACDM 
cosmology with the following model parameters: z r = 20, t* 
3 Gyr, a, = 0, T w = 2 x 10'' 
were computed using Mo 
corresponds to the universe as a whole 



6 K and rfw = 1.5. These fractions 
+00 and zq — —1, and therefore 



cold gas is ejected from the disc with a typical wind ve- 
locity u w and with a typical outflow rate 

(26) 



~ 1 — 5, the wind 
L the wind veloc- 



The two additional parameters are rj^ 
effi ciency, and m w ~ 200 — 500 km s" 
ity l|Springel fc HernauistlE003rj|) . These wind parameters 
are typical of observed outflows in star forming galaxies 
llMartmlll999l . 

The fate of this ejected gas depends on the halo mass. 
If the wind velocity exceeds the escape velocity of the halo, 
the ejected gas leaves the halo into the diffuse background, 
from where, eventually, it will be accreted again. Such 
winds are referred to as 'unbound'. If, on the other hand, 
the halo is too massive, the ejected gas remains in the hot 
halo component, from where it will eventually cool again. 
Such winds are referred to as 'bound'. 

Assuming again that all baryons are locked up into the 
cold component, we can compute the global wind outflow 
rate by integrating over the EPS mass function. We finally 
obtain the following equations, valid in the general case 

/wind = ^w/cold- (27) 

where the outflow rate, in units of Gyr -1 , is given by 

= i)wW,. (28) 

We first compute the unbound fraction. It corresponds to 
winds emitted by halos whose escape velocity is smaller 
than the wind velocity. Using the EPS distribution, we 
get 

erfc(t' min /v / 2) - erfc(i/ w /\/2) 



where v w is defined by the 'Wind Mass' 
5 c (t) - S c (t Q ) 



v w (M ,to) 



V^(M W ) 2 - ct(M ) 2 



(30) 



erfc(i/ min /v / 2) 



(29) 



This Wind Mass' (the halo mass above which winds are 
bound) is related to the wind velocity by noticing that for 
typical dark matter halos, u CS c — 3V2oo- Using the stan- 
dard Virial relation (Eq. EJ), we obtain the wind Virial 
temperature 

feTw = -^FiH^. (31) 

lo 

The bound fraction is just 1 — £ w . The fate of the un- 
bound gas depends now on the parent halo mass M . If 
M > M w , the gas is recycled into the halo diffuse com- 
ponent, and ultimately into the halo hot component as 
z — > zq. If Mq < M w , the gas is lost into the intergalactic 
medium, outside the boundaries of the parent halo, and 
never come back. Note that in the latter case, £ w is al- 
ways equal to one. This very crude model turns out to 
be surprisingly accurate in predicting results obtained by 
numerical simulations (see the following sections). 



3. 7. The Baryon Supply Chain 

We are now in a position to compute the baryon bud- 
get history. The last sections were devoted to computing 
mass transfer rates between our 4 baryon components. 
The mass fraction in each component are the independent 
variables in our problem: / bac k, hot, /cold and /», referring 
respectively to diffuse background, hot gas, cold discs and 
stars. 
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reionization redshift 


U 


star formation time scale 


»7w 


wind efficiency 


Cw 


unbound wind fraction 


/. 


stellar fraction 


/cold 


cold gas fraction 


/hot 


hot gas fraction 


/back 


background gas fraction 




star formation rate 


k-'cool 


cooling rate 


^acc 


accretion rate 



Table 4. Main notations 



The methodology we follow in this paper allows us 
to compute in advance 3 important mass transfer rates. 
This rates are w acc , u> coo i and w*, referring respectively 
to the Cosmic Accretion Rate, the Cosmic Cooling Rate 
and the Global Star Formation Rate. Our very crude wind 
model provides us with an additional parameter, namely 
the unbound fraction £ w . These various rates are plotted 
in Figure for our fiducial model and our notations are 
summarized in table l3~TI 

We have to solve a set of ordinary differential equa- 
tions, with pre-computed transition rates between each 
component of our baryon supply chain. 

— "77 — = Cw^w^/cold — 'J-'acc/back, (32) 



^acc/back ^cool /hot + (1 - Cw)?7w^*/coid, (33) 



4. Simulations versus Model 

We now compare our analytical predictions to the baryon 
budget history obtained in our high-resolution hydrody- 
namical simulations. Recall that we can compute analyt- 
ically the baryon history for the universe as a whole, but 
also on a halo by halo basis. In order to make a care- 
ful comparison, we need to extract halos from the simu- 
lated density field. We use for that purpose a halo detec- 
tion code based on the Spherical Overdensity algorithm 
l Lacev fc Coldllfffl^ . We also need to carefully define the 
effective mass resolution of our simulations, with respect 
to star formation. This additional mass scale is a pure nu- 
merical artifact, that can be accounted for explicitly in 
our analytical model. Using this modified model, we will 
estimate how our results converge (or not) to the correct 
halo model predictions. 

We need also to estimate the halo star formation time, 
(as defined in the previous section) in our numerical simu- 
lations. The star formation algorithm is based on a simple 
Schmidt law, with a specified (over-) density threshold. It 
is however more complex than the approach used in the 
halo model. We will show that both approaches can be 
related to each other, with a 'shape factor' reflecting the 
probability distribution function of the cold phase density. 

We finally need to estimate the unique free parameter 
in our analytical model: the orbital decay timescale. This 
sort of 'model calibration' will be performed directly using 
our simulation results. We will also extend this 'model ver- 
sus simulation' c omparison to other numerical d ata kindly 
provided to us bv lSrjringel fc Hernau ist ( 2003bJ), for which 
galactic winds were included. 



cold 



dt 



^cool /hot — ^*/c. 



?7wW*/cold, 



df* _ , 

—7T — J cold- 

dt 



(34) 



(35) 



If Mo > M w , one sees that the total baryon mass is con- 
served. Note that Eg. 1321 should be modified if Mq < M w : 
the wind contribution is set to and therefore the total 
baryon mass in the parent halo is not conserved anymore. 
This is as expected, since winds are now escaping outside 
the parent halo boundaries. 

Using any time integration method of sufficient accu- 
racy, one can finally solve for the previous set of differential 
equations. We used in this work a Backward Euler scheme. 
Interestingly, one can solve formally the latter system us- 
ing matrix exponentials. This type of equations ar e typical 
of gala xy formation stu dies, like the ear ly work of lTinslevI 
l)l980|) . More recently, IPei et alJ {l999) have designed a 
similar approach, based on the observed galaxy luminos- 
ity functions, while here, our equations are based on EPS 
theory. Figure shows the baryon budget evolution for 
our fiducial model. Before applying this analytical model 
to cosmological observations, we need to determine its va- 
lidity range using high resolution numerical simulations. 



4.1. Halo Detection 

It is an absolute necessity to define a halo in a numer- 
ical simulation as it is defined in the theoretical model. 
As explained before, in order for the Press & Schechter 
approach to be valid, the halo mass is the mass enclosed 
in radius i?200i enclosing an overdensity 200 times larger 
than the av erage background density. As no ted by sev- 
eral authors l|jenkins et al.ll200ltlwhitell2002() . this rather 
large region encloses dark matter particles which are not 
completely relaxed yet, and also several large satellites fly- 
ing by. As a consequence, the halo mass turns out to be 
highly dependent on the exact algorithm used to detect 
automatically dark matter halos in the simulation. 

These same authors suggest the following strategy: in 
a first step, halos are detected with a high density contrast 
(A = 600 is our choice here) by any classical algorithm (we 
use Spherical Overdensity in this paper) . Since the density 
contrast is high, the detected region is in a well relaxed 
state and all algorithms agree more or less on the mass 
and on the number of detected halos. The halo center is 
defined as the center of mass of this high density region 
only. In a second step, the halo radius is increased up to 
-R200, in order to obtain the large halo mass required by 
the Press & Schechter prediction. 
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We then compute for each individual halo the total 
stellar mass within i?200- Since we have stored the for- 
mation epoch of each individual star particle, we can also 
compute the complete star formation history of the parent 
halo. This last point is very important: we do not compute 
star formation rates and gas content of individual galax- 
ies. Since each halo can host several galaxies (one central 
and several satellites), the galaxy baryon budget will be 
somewhat different than the halo baryon budget. Within 
the halo radius, we also compute the fraction of cold gas, 
defined as T < 10 5 K, and p > 10 5 p. The remaining gas is 
considered as 'hot gas', even though its temperature can 
be lower than the Virial temperature of the halo. 

The effective star formation rate of each halo is com- 
puted by simply dividing the total amount of star created 
during the last 10% of the halo age by the elapsed time. 
The results presented in this section will be based on this 
analysis. Each individual halo baryon budget are averaged 
into mass bins, in order to compare with the halo model 
predictions. 

4.2. Mass Resolution 

Results of cosmological simulations depend strongly on 
the mass resolution of the code. The two main numerical 
limitations are the box length and the number of particles. 
As shown in Figure^U the more particles we start with, the 
more small mass halos and galaxy satellites we obtain in 
the simulations. Since we are interested here in the global 
baryon budget, each individual halo must be able to allow 
gas to cool down and condense, and ultimately form stars. 
This is a more stringent requirement than just reach the 
Virial overdensity of A = 200. 

Since star formation occurs at the high end of the 
density distribution, we take the star formation density 
threshold as the limiting factor that defines our mass res- 
olution. Let us assume for sake of simplicity that each 
halo is a pure isothermal sphere. The gas density profile 
is given by 

-2 



0.10 - 



-p 



R 



2011 



(36) 



10 



7/V, 



(37) 



The radius above which star formation occurs is given by 
p(r) > po(l + z) a ° . If we require that within this radius, 
we have at least 10 dark matter particles of mass m p , we 
obtain the minimal halo mass as 

p (l + z) Q ° 1 1/2 
A/3p(z) _ 

One clearly sees that the higher the density threshold 
for star formation, the larger the minimal mass will be. 
For star formation density thresholds constant in co- 
moving units (ao = 3), this mass scale is a constant 
in time. This is the case for the AMR simulations pre- 
sented here, for which, using Equation 1371 we obtain 
M roso i ~ 400m p . On the other hand, for star formation 
density threshold constant in physical units, the mass res- 
olution scales as (1 + z) -15 . SPH simulations presented 
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Fig. 8. Time evolution of the various baryon phases in our 
highest resolution run L10N512S30. Symbols are mass frac- 
tion measured in the simulation (squares: diffuse background, 
crosses: hot gas, diamonds: cold gas and triangles: stars). Lines 
are the predictions of our analytical model with A/ roso i ~ 2x 10 s 
h -1 M0 (solid: diffuse background, dot-dot-dashed: hot gas, 
dot-dashed: cold gas and dashed: stars). 



in ISpringel k, HernauistJ l)2003bj) were based on this sec- 
ond approach. The mass resolution we obtain in this 
case is M reso i ~ 1000(1 + z) -1 - 5 m p . At high redshift, 
z ~ 20, the corresponding mass resolution can be as low 
as M reso i ~ 10m p . 

Simulated halos with mass lower than -Mresoi will not 
be able to form stars or, equivalently, condensed cool gas. 
Therefore, they will be part of the simulated diffuse back- 
ground. This new mass scale is a pure numerical arti- 
fact, that strongly affects our results. We take this mass 
scale into account in our analytical model by setting the 
Minimal Mass for star forming halos M m ; n as the maxi- 
mum between the true physical Minimal Mass and M rcso i. 
As we will see later in this section, this trick is a very 
powerful tool to account for finite resolution effect in the 
simulation, and to assess the convergence properties of our 
numerical results. 

4.3. Halo Star Formation Time 

The methodology we use in this paper is to describe star 
formation on a halo by halo basis. We completely discard 
the detailed modeling of exponential gas discs and nuclear 
bursts. This is usually performed in semi-analytical mod- 
els of galaxy formation. In our AMR simulations, we do 
have however a higher level of complexity than in the an- 
alytical model. Visual inspection of density maps shows 
the presence of gas discs in centrifugal equilibrium, as 
well as several small satellites orbiting a central galaxy 
(see Fig. nj. We are aware of the fact that many physical 
ingredient are probably missing in our current numeri- 
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Fig. 9. Average star formation rate of simulated halos in var- 
ious Virial temperature bins, in unit of M co id/t*(pt)- in om " 
framework, this is a direct measure of the 'shape factor' F(p) 
of the underlying cold gas density distribution. Numerical data 
suggest a constant value represented here as the dashed line 
F(p) ~ 3. Diamonds are for run L1N256S30 at z = 5.5, tri- 
angles for run L10N256S30 at z = 2.5 and squares for run 
L100N256S30 at z = 0. 
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Fig. 10. Average star formation rate of simulated halos in 
various Virial temperature bins, in unit of Mh ot -R200/V200 for 
the high efficiency runs. This is a direct measure of the cool- 
ing rate of hot gas into dense cold discs. At high tempera- 
ture (T > 10 7 K), we observe a sudden drop due to inefficient 
Bremstrhalung cooling. At lower temperatures, the cooling rate 
has a plateau around 1/3. In our framework, this suggests 
an orbital decay timescale i or b — 3.R200/V200 for infalling gas 
clumps. Diamonds are for run L1N256S3 at 2 = 5.5, triangles 
for run L10N256S3 at z = 2.5 and squares for run L100N256S3 
at z — 2.5 also. 



cal solution of galaxy formation. Nevertheless, we need to 
establish the link between our analytical model and our 
numerical implementation of star formation. 

Since star formation in the code is based on a Schmidt 
law (see Eq. QJ, we can compute the instantaneous star 
formation rate in any halo by integrating over the entire 
cold gas present in that particular halo. 



M , = M, 



cold 



t*(p) 



(38) 



where pt = p${l-\-z) a ° is the star formation density thresh- 
old and p{p) is the mass fraction of cold gas with density 
p. The exact form of the cold gas distribution function p 
is beyond the scope of this paper. It is likely to be deter- 
mined by the global surface density as well as the small 
scale turbulence inside rotating discs. We will make here 
the very crude approximation that p is self-similar in the 
variable p/pt, so we can simplify the last equation further 
more into 



M* = 



Meld 

t*(pt) 



F(P), 



(39) 



where F (p) is a dimensionless 'shape factor' that depends 
on the exact form of the cold gas density distribution 



/oo 
p{x)x 1,2 dx. 



(40) 



We are now in a position to make a direct link between our 
analytical model and numerical simulations. We recognize 
in the last equation the halo star formation rate as defined 
in section 13.51 with star formation parameters given by 
a* = «o and t* = to/F(p). The shape factor F has the 
effect of reducing the effective halo star formation time 
scale, relative to the reference time to. Indeed, if very high 
density gas is present, the star formation rate is likely to 
increase accordingly. 

Since we can't predict the value of this shape factor, 
we have to measure it directly in the simulations. We 
plot in Figure [5] the halo star formation rate, in units of 
Af co id/£*(/9t) ( see Eq.Oni), and averaged over halos of sim- 
ilar mass. This should be equal to the shape factor F(p). 
For 3 different box sizes and at 3 different redshift, this fac- 
tor is not exactly a constant, although it varies slowly with 
mass. This illustrates that our approach is only a first or- 
der approximation of our simulation results. Nevertheless, 
we approximate this by taking F(p) ~ 3, as suggested by 
the dashed line in FigureEl This specifies how star forma- 
tion in the simulations and star formation in the model 
are connected to each other. 

4.4. Halo orbital decay timescale 

The only unknown parameter in our analytical model is 
the orbital decay timescale of infalling gas clumps (see 
appendix [UJ , before they reach the high-density disc in 
the halo center. When cooling is very fast, for halos with 
Virial temperature T min < T200 < T max , we have assumed 
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Fig. 11. Global comoving star formation rate as a function of redshift. In each plot, symbols are for numerical simulations, while 
lines are for our corresponding analytical prediction. Upper left plot: high efficiency series and low efficiency series with run 
L1N256S3 (black squares), run L1N256S30 (grey squares), run L10N256S3 (black triangles), run L10N256S30 (grey triangles), 
run L100N256S3 (black diamonds) and run L100N256S30 (grey diamonds) . Lower left plot: convergence study series with run 
L10 N512S30 (diamonds), run L10N 256S30 (triangles), run L10N128S30 (squares) and run L10N64S30 (crosses). Upper right 
plot: \Svrinael & Hern ouisi Ii 2003&) simulation s with run 03 without winds (triangles) and Q3 including winds (diamonds) . 
Lower right vlot ^Svrkiael & Hernauisl \200SA) simulations including winds with run Q5 (stars), run Q4 (diamonds), run Q3 
(triangles), run Q2 (squares) and run Ql (crosses). The solid line is the 'fully converged' model prediction. 



that the accretion rate into the disc is controlled by the 
orbital time scale of infalling satellites. Computing this 
time scale is beyond the scope of this paper. It is prob- 
ably determined by details in the gravitational dynamics 
and satellite dynamical friction. These aspects are all key 
ingredients of semi-analytical models of galaxy formation. 

In order to determine this orbital decay timescale, we 
perform again a direct analysis of our numerical simula- 
tions. Let us consider the case of very fast star formation 



to = 3 Gyr and ovo = 3. In this case, the halo star forma- 
tion rate is almost equal (within 10%) to the halo cooling 
rate. This can be later confirmed by the analytical model. 
We plot in Figure[!|]the average star formation rate of ha- 
los within different mass range, in units of A/hot ^200/^200- 
In our framework, this quantity is exactly equal to the 
ratio (i?20o/V20o)/^orb- Here again, this ratio is not per- 
fectly a constant, illustrating the fact that our model is 
only a first order approximation, but for the 3 different 
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box sizes and at 3 different redshifts, the curve exhibits a 
plateau around i or b — 3i?2oo/^20o- We take this value as 
our canonical value in the analytical model. 

4.5. Global Baryon Budget 

We now present in greater details our simulation results, 
starting with the baryon history for the universe as a 
whole. Figure |S] shows the baryon history in our high- 
est resolution run L10N512S30. The run parameters cor- 
respond to a low efficiency star formation model. Each 
phase is defined by well defined limits in the p - T dia- 
gram, as defined in Section 12.31 The various symbols in 
Figure |S| refer to baryon fractions in different snapshots 
of the simulation, while lines refer to the analytical model 
predictions, with M roso i ~ 2 x 10 s h _1 M Q , as given by 
Equation E| The other parameters of the model are set 
to their standard values {F(p) = 3 and i? or b = 3i?2oo)- 
The agreement between the simulation and the model is 
very good (within a factor of 2), given the simplicity of 
the latter. 

We now examine more closely the global star formation 
rate as a function of redshift, measured in all our simu- 
lations, and compare our various results to the analytical 
model. This quantity is a key prediction of hierarchical 
model of galaxy formation. It translates more or less di- 
rectly into galaxy colors and luminosities, and provides a 
stringent test of the current cosmological theory. As star 
particles are created during the course of a simulation, we 
keep track of their birth epoch. It is straightforward to 
compute, using the last output only, the star formation 
epoch histogram (history). 

Figure ^] shows this global star formation history 
for our 'convergence study' simulation suite: L10N64S30, 
L10N128S30, L10N256S30 and L10N512S30. Numerical 
results are shown as symbols, while analytical predictions 
are shown as lines. The analytical model predictions are 
computed with a Minimal Mass corresponding to the mass 
resolution of each run, as given by Eauationl37l Since our 
star formation recipe is based on a constant overdensity 
threshold, the mass resolution is AT rcso i ~ 400m p . The 
solid line stands for the analytical prediction, without any 
finite resolution effects (Af roso i ~ 0). This gives an indica- 
tion on how our results have converged to the 'true' star 
formation history in this particular model. 

Around z ~ 3 — 5, our highest resolution run 
L10N512S30 is close to the correct value. At very high 
redshift however, the star formation rate is lower than the 
expected value by a factor of ten. The mass resolution 
-Mrcsoi is indeed significantly higher than M m i n at redshift 
z > 10: this explains the origin of this discrepancy. As il- 
lustrated by the Figure EH low resolution runs do miss the 
formation of dwarf galaxies that contribute significantly 
to the global star formation history. 

This first series of simulations was performed for a box 
length L = 10 h _1 Mpc. By z ~ 3, the non-linear scale 
becomes comparable to the box size. We have to stop the 



simulations, as our realizations are not representative of 
the universe as a whole anymore. 

Our second series of simulations was designed to ex- 
plore larger (L — 100 h^ 1 Mpc) and smaller (L = 1 h _1 
Mpc) scales, as well as another, more efficient, star forma- 
tion scenario with to — 3 Gyr. The corresponding Schmidt 
law is therefore ten times more effective than in the first 
case. The resolution was fixed to 256 3 particles. All six re- 
sults are shown in Figure lTTl with symbols, while the corre- 
sponding analytical predictions are overploted with lines. 
Note that the smaller box (I = 1 h _1 Mpc) has a mass 
resolution smaller than M m j n . It has therefore converged 
to the 'true' star formation history. The analytical model 
is indeed in good agreement with our numerical results. 
This very small scale simulations have to be stopped even 
earlier than the previous ones (z ~ 5). The largest box 
size (L = 100 h _1 Mpc), on the other hand, is strongly af- 
fected by finite mass resolution effect. The star formation 
history is drastically different than the other two. Star for- 
mation begins very late, around z ~ 5, and the peak value 
is one order of magnitude lower than the 'true' star for- 
mation rate. Our analytical model provides again a good 
fit to numerical data, when the poor mass resolution is 
taken into account to define the diffuse background. 

Let us now compare the two different star formation 
scenario (inefficient with to = 30 Gyr and efficient with 
to = 3 Gyr). Both star formation history are parallel to 
each other, but with a factor of two difference only. For 
the very efficient scenario, star formation is limited by ac- 
cretion and cooling, rather than by the Schmidt law. As 
discussed in the next section, very efficient star formation 
models give SFR quasi independent of the physical param- 
eters. This is confirmed by our numerical simulations. The 
agreement between the simulation and the model is ex- 
tremely good, except at low redshift (z < 3), for our large 
box runs, where the model significantly overestimates the 
star formation rate. This disagreement might be reduced 
by improving the analytical model, along several lines we 
have outlined in this paper. Note that finite volume ef- 
fect might also have an additional effect on the numerical 
simulation predictions, but we do not try to include those 
subtleties in the analytical model. 

From now on, we have analyzed AMR simulations 
for which star formation is triggered in region where the 
gas density exceeds an overdensity threshold. In this case, 
we naturally compare our numerical results to the 'ac- 
celerated efficiency' analytical model. We now test the 
model predictions for the 'constant efficiency' analyti- 
cal model, using SPH result s kindly provided to us by 
ISnringel fc HernciuistJ l|2003bh . Our interest to these SPH 
results is twofold: first, star formation is triggered in 
regions where the gas density exceeds a physical den- 
sity threshold. This scenario, as already discussed in 
Section ^31 correspon ds to a 'constant efficiency' sta r for- 
mation model. Second. lSnringel fc Hernauistl l|2003b|) have 
included galactic winds in their numerical modeling, giv- 
ing us the opportunity to test our simple feedback model. 
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Fig. 12. Halo star formation rate, in units of A'hoo/t, as a function of the halo Virial temperature. In each plot, symbols are 
for numerical simulations, while lines are for our corresponding analytical prediction. Upper left plot: high efficiency series with 
run L1N256S3 at z = 5.5 (diamonds), run L10N256S3 at z = 2.5 (triangles) and run L100N256S3 at Z = 0. Lower left plot- 
low efficiency series with run L1N256S30 at z — 5.5 (diamonds), run L10N256S30 at z = 2.5 (triangles) and run L100N256S30 
at z — 0. Lower right plot: convergence study series at z = 2.8 with L10N64S30 (diamonds), run L10N128S30 (triangles), run 
L10N256S 30 (squares ) and L 10N512S30 (crosses). The solid line is the 'fully converged' model prediction. Upper right plot: 
\Svrinael & Hernauisi \200SA) simulations including winds, with run R4 (diamonds) at z = 6, run Q4 at z = 3 (triangles) and 
run D4 at 2 = 1.5 (squares). 



In Figure ^] are shown ISprineel h HernauistJ (|2003b^ 
results when galactic winds are turned off. The only dif- 
ference with our simulations, apart from the overall nu- 
merical techniques, comes from the star formation details. 
They have considered the following parameters to = 2.1 
Gyr, iiq = 0.1 cm~ 3 and «o = 0. For the analytical model, 
we use the same 'shape factor' F(f£) — 3, as determined 
earlier using our AMR results. This translates into the fol- 
lowing parameters for the analytical model: t* = 0.7 Gyr 



and a* = 0. Taking into account the finite mass resolution 
(using Equation II Q , we can now compare our analytical 
prediction to SPH results. The agreement is clearly very 
good: this is rather encouraging for our model, since our 
main unknown parameter, the mean orbital length, was 
kept fixed to its canonical value i? or b = 3i?2oo , calibrated 
on our AMR results. 

We now analyze SPH results when galactic winds are 
turned on. Recall that winds are assumed to eject cold 
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Fig. 13. Mass fraction in each baryon phase as a function of halo Virial temperature for the high efficiency simulation series. 
In each plot, diamonds refer to run L1N256S3 at z = 5.5, triangles to run L10N256S3 at z = 2.5 and squares to run L100N256S3 
at * = 0. Lines are for the corresponding analytical model. The upper left plot shows the star mass fraction; the lower left plot 
shows the hot gas mass fraction; the lower right plot shows the total baryon fraction and the upper right plot shows the cold 
gas fraction. 



gas from the rotating discs at a typical wind velocity u w ~ 
250—500 km/s and with an efficienc y parameterized by 77 w . 
This ap proach, directly inspired bv lSpringel fc Hernauistl 
(2003b), allows a straightforward comparison between our 
analytical model and SPH results with winds. This com- 
parison is shown in Figure ^2 The analytical model pa- 
rameters were set to T w — 2 x 10 6 K and rj w = 3. The 
Virial temperature T w (below which winds are unbound 
from their parent halo) correspon ds closely to the wind 
velocity u w ~ 500 km/s used by ISorineel fc Hernauistl 
l|2003blj) in their SPH simulation (see Ea.l3T)l, The wind ef- 
ficiency parameter, however, was chosen 50% higher than 



the value 77 w = 2 used in the SPH sim ulation. As sug- 
gested bv lSpringel k. Hernauistl |2003b), the global wind 
efficiency (at the halo scale) is higher than the local wind 
efficiency (at the star forming regions scale) due to gas en- 
trainement: additional cold gas, lying outside star forming 
regions from which winds originate, can be expelled by the 
ram pressure of the outflow. Figure ITU shows again a good 
agreement between SPH simulation results and our ana- 
lytical model. 

We consider now SPH re sults from the 'Q series' in 
ISorineel fc Hernauistl ll2003bl) paper, in order to study fi- 
nite mass resolution effects. Wind parameters are fixed 
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Fig. 14. Mass fraction in each baryon phase as a function of halo Virial temperature for the low efficiency simulation series. In 
each plot, diamonds refer to run L1N256S30 at z = 5.5, triangles to run L10N256S30 at z = 2.5 and squares to run L100N256S30 
at z = 0. Lines are for the corresponding analytical model. The upper left plot shows the star mass fraction; the lower left plot 
shows the hot gas mass fraction; the lower right plot shows the total baryon fraction and the upper right plot shows the cold 
gas fraction. 



to T w = 2 x 10 6 K and rj w = 3. Note that in this 
case, due to a star formation strategy based on a con- 
stant physical density, M roso i is now varying with time 
M rfiso i ~ 1000(1 + zy-^Wp . In Figure El are plotted 
ISpringel fc Herncmistl l|2003rjj) simulation results together 
with our model predictions, when mass resolution effects 
are taken into account. At high redshift, we recover the 
correct convergence towards the asymptotic (M roS oi = 0) 
converged curve. At intermediate redshift, however, the 
agreement is getting worse, although both curves remain 
close to each other within a factor of 2. One possible rea- 
son is that Equation|^|is not accurate enough to estimate 



the effective halo mass resolution of SPH simulations, es- 
pecially in presence of winds. 

We conclude that the global baryon history we ob- 
tain in numerical simulations is correctly described by our 
analytical model, when finite mass resolution effects are 
taken into account. This was tested for 2 different nu- 
merical technics (AMR and SPH), two different star for- 
mation scenarios ('constant' and 'accelerated' efficiency), 
with and without galactic winds. The main effect of insuffi- 
cient mass resolution is to artificially decrease the average 
age of stars in the universe and to lower the star formation 
rate at peak value. 
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Fig. 15. Mass fraction in each baryon phase as a function of halo Virial temperature for the convergence study simulation series 
at z = 2.8. In each plot, diamonds refer to run L10N64S30, triangles to to run L10N128S30, squares to run L10N256S30 and 
crosses to run L10N512S30. Lines are for the corresponding analytical model. The upper left plot shows the star mass fraction; 
the lower left plot shows the hot gas mass fraction; the lower right plot shows the total baryon fraction and the upper right plot 
shows the cold gas fraction. The solid lines are the 'fully converged' model predictions. 



4.6. Halo Baryon Budget 

We present now our simulation results concerning the 
baryon budget in individual halos. Halos are detected 
in the simulations according to the method explained in 
Section O The Extended Press & Schechter (EPS) the- 
ory give us the opportunity to apply our analytical method 
to an 'average' halo of mass Mq at redshift zq. The model 
predictions can then be compared to the average baryon 
fractions (in each of the 4 different baryon phases), the 
average being taken over all halos in a given mass range, 
centered around Mq. 



Using all star particles found within i?200, we com- 
pute the star formation history within each halo. We 
then compute the star formation rate as explained in 
Section 14.11 Figure 1121 shows our various numerical re- 
sults, in cluding the SPH simulations w ith winds provided 
to us bv lSpringel fc Hernauistl l)2003bj) . together with the 
analytical model predictions. We have plotted the halo 
star formation rate as a function of the halo tempera- 
ture T200, m un it s °f Mvaa/t. This definition corresponds 
closely to the 'specif i c star formation rate' defined by 
ISprineel fc Hernauistl l|2003bj) . 
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Each curve shows a sharp cut-off at the low mass end, 
corresponding to the Minimal Mass for each run. The 'con- 
vergence study' series clearly illustrates that this Minimal 
Mass is in fact equal to the mass resolution M reso i, ex- 
cept for our smallest box size. Another non trivial effect 
of finite mass resolution is to increase the halo star forma- 
tion rate. This is due to a higher Cosmic Accretion Rate 
(/acc oc cr(M m i n ) _1 , see Eq. IT5)l . as M m i n is increased. 
Our analytical model predictions are in good agreement 
with the halo star formation rates, even when winds are 
present, as soon as finite mass resolution is explicitly ac- 
counted for in the model. The role of winds is to remove 
cold gas from small mass halos T200 < 2 x 10 6 K, so that 
the halo star formation rate decrease accordingly. 

When we compare the high star formation efficiency 
series with t = 3 Gyr with the low star formation effi- 
ciency series with to = 30 Gyr, we see that the former 
has a higher halo star formation rate at high redshift, but 
a much lower halo star formation rate at low redshift, as 
cold gas is almost completely consumed. This complex be- 
havior is well reproduced by the analytical model. Both 
series show a sharp decline of the halo star formation rate 
at the high mass end: the cooling efficiency decreases for 
high mass halos and the mass accretion rate on cold discs 
vanishes. Note that we assumed here a zero metallicity 
plasma. Gas cooling is likely to be more efficient around 
T200 — 10 7 K if metals are present. 

We now present the baryon budget inside dark matter 
halos, as a function of the halo Virial temperature. In each 
halo, we compute the total baryon fraction, which can be 
further decomposed into hot gas, cold gas and stars, using 
the definitions of Section l2.3l and l4.ll Figure [R3l shows this 
halo baryon budget for the 'high efficiency' simulation se- 
ries, FigurelTHfor the 'low efficiency' simulation series and 
Figure IT51 for the 'convergence study' simulation series. 

In each case, there is good agreement between numer- 
ical and analytical results. The total baryon fraction is 
close to the average value /b for halos more massive than 



1.000 r 



M„ 



and vanishes for smaller halos. Note that the to- 



tal baryon fraction is actually slightly above the universal 
value for massive halos. Dissipative collapse of baryons 
condense more mass than collisionless collapse of dark 
matter. The analytical model predicts a sharp transition 
for M < M reso i where the total baryon fraction vanishes. 
In the simulations, this transition is much smoother, but 
its location is correctly predicted around M reso i. On the 
other hand, when the mass resolution is small enough, the 
low mass end of the analytical curve has also a smooth 
transition. This comes from halos whose mass was greater 
than M m ; n at early time, but as reionization heats up the 
background temperature, these same halos end up with a 
mass smaller than Af m i n . The total baryon fraction at the 
final redshift is therefore smaller than the universal value 
but greater than zero. 

The hot gas fraction is also strongly affected by reso- 
lution effect. As we resolve smaller and smaller mass pro- 
genitors, more and more gas is condensed into cold discs 
and eventually stars. Figure ED shows the hot gas fraction 
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Fig. 16. Cosmic Star Formation Rate for different model pa- 
rameters: t, = and n v = 0, z T = 20 (solid line), z r = 12 
(dotted line), z r = 6 (dashed line). 



for the 'convergence study' series and one clearly sees that 
our poorest resolution run overestimates the hot gas frac- 
tion by a factor of 3. As a consequence, the cold gas and 
star mass fraction gradually increase as the mass resolu- 
tion is improved. The curve showing the cold gas fraction 
as a function of Virial temperature is very similar to the 
halo star formation rate, as expected from our analytical 
model, although they have been computed completely dif- 
ferently in the simulation analysis. In the 'high efficiency' 
series, our largest box run L100N256S3 has almost en- 
tirely consumed the cold gas of the most massive halos, 
where cooling has virtually stopped (see Fig. 113)1 . The cor- 
responding run with a 'low efficiency' star formation sce- 
nario still shows some surviving cold gas at the high mass 
end of the halo distribution. 

We conclude at the end of this section that our an- 
alytical model has proven to successfully reproduce the 
complex b ehavior observed in our simulations and in sim- 
ulations of ISpringel fc Hernauistl l)2003bj) (with winds) if 
finite resolution effects are taken into account. We note 
however a slight tendency of the model to overestimate 
the star formation rate at low redshift. The baryon bud- 
get was analyzed in great details on a halo by halo basis 
(and in an average sense). We emphasized the important 
role played by the Minimal Mass M mia in controlling the 
baryon history in each mass range. We also described how 
cooling and winds introduce characteristic features in the 
various plots showing the baryon evolution as a function 
of the halo Virial temperature. 

5. Observations versus Model 

From now on, we assume that our analytical model gives 
an accurate (within a factor of 2) modeling of the baryon 
history in a hierarchical universe. The careful comparison 
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Fig. 17. Cosmic Star Formation Rate for different model parameters. Left plot: z r = 20 and ?y w = 0, i* = (solid line), t* = 0.1 
Gyr (dotted line), t* = 1 Gyr (dashed line), t* — 3 Gyr (dot-dashed line), t» = 10 Gyr (dot-dot-dashed line). Right plot: 
z r = 20 and U = 3 Gyr, rj w = (solid line), r}„ =0.1 (dotted line), n v = 1 (dashed line), r)„ = 3 (dot-dashed line), rj„ = 10 
(dot-dot-dashed line). 



to numerical simulations we performed in the previous sec- 
tion was a necessary step to calibrate our model and to 
estimate its level of accuracy. We are now in a position to 
use this model as a tool to analyze current observational 
constraints put on the baryon history in the universe. 

5.1. Model parameters study 

We briefly recall how the analytical predictions depend on 
the model parameters. We consider that all cosmological 
parameters are fixed to their 'concordance' ACDM value, 
as we did throughout this paper. We also assume that the 
reionization temperature is given by T r ~ 6 x 10 3 K. We 
finally use the cooling model of Section [UJ valid for a pri- 
mordial, zero metallicity, H and He plasma and a wind 
velocity u w ~ 500 km/s. We end up with 3 main param- 
eters that we allow to vary in our model: the reionization 
redshift z T , the star formation time scale (or consump- 
tion time scale) i* = M co id/M* and the wind efficiency 
n w = M w i n( i/M*. We consider in this section only 'con- 
stant efficiency' star formation models a* = 0. The alter- 
nate scenario with 'accelerated efficiency' star formation 
and a* = 3 gives similar, though slightly poorer, results. 
Using current observational constraints, it is not easy to 
discriminate between these two models. 

We now present the model predictions for the Cosmic 
Star Formation Rate, for different values of z T , i* and rj w . 
Figure 1161 corresponds to an 'infinite efficiency' star for- 
mation model (f» = Gyr), for various reionization red- 
shifts. In this extreme case, the cosmic star formation rate 
(CSFR) is equal to the mass accretion rate of cold discs, 
since star formation is instantaneous. At reionization, the 
background gas is promptly heated to 10 4 K and the star 



formation rate drops. We notice that at redshifts lower 
than 6, all models agree with one another, leading us to 
the conclusion that, as soon as reionization proceeds early 
enough, this parameter is in fact irrelevant to the low red- 
shift universe. 

Figure 1171 shows the influence of the two remaining 
(and really interesting) parameters. Star formation inside 
cold, centrifugally supported, galactic discs acts merely as 
a delay with respect to diffuse gas accretion. The epoch 
of peak star formation is delayed from z ~ 4 for t* = 
down to z ~ 1 for £„ = 10 Gyr. It is worth noticing that 
the amplitude of the global star formation rate is mainly 
determined by the cosmological accretion rate, and that 
the small scale, poorly known, physics of star formation 
introduces a rather small modification to the curve. When 
winds are included in the model, they lower significantly 
the amplitude of the CSFR. They also have the non trivial 
effect of advancing the epoch of peak star formation, from 
z ~ 1 for r) w = up to z ~ 3 for r/ w = 10. They play an 
important role at low redshift whereas at high redshift the 
CSFR seems independent of ry w . Indeed, for the 'constant 
efficiency' recipe the star formation rate at high redshift is 
so small compared to the accretion rate that winds (which 
are proportional to the star formation rate) doesn't affect 
the amount of cold gas. 

Computing the baryon history using our analytical 
model is a very fast operation. This gives us the opportu- 
nity to compute various observational quantities for a grid 
of model parameters. Using various observations, we will 
now try to constrain our parameter space and shed light 
to these 2 important galaxy formation ingredients: 4* and 
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5.2. Cosmic Star Formation Rate 

One of the main goal of this paper is to compute 
the star formation history in a hierarchical universe. 
We will now compare the model predictions to the 
observed star formation rate. Figure ^] shows the 
observational data points, usually referr ed to as the 
'Madau plot', compiled by lElbad l)2005l) and unifor- 
mally corrected from cosmological distances, incom- 
pletness an d dust abso r ption Original data point s 
came from iHughes et all lll99Sl): ISteidel et all JlQQflh: 
iFlores et al.l il999h: lOlazebrook et al.l Jl999h: lYa.r, et al.l 
l|l999|k iMassarotti et all l|200l|) : iGiavalisco et alJ (|200J). 
These data points indicates an epoch of peak star forma- 
tion rate at z ~ 2, followed by a rapid fall off (by a factor 
around ten) between z = 1 and z = 0. 

The shape and the normalization of the Madau plot 
are well reproduced by the parameter choice ~ 1.5 
Gyr and r/ w ~ 1. It is worth noticing that this star for- 
mation time scale corresponds roughly to the value one 
can infer from the consumption timescale in local galaxies 
l|Kennicutt et alJ Il994t iKennicuttl Il998|) . The wind effi- 
ciency we obtain from this fitting exercise is close to the 
value inferred from H a observations of high-SFR galaxies 
ljMartirJll999I) . Galactic winds are required in order to re- 
produce the rapid fall off of the star formation rate at low 
redshift. 

The observational constraints on our two main star 
formation parameters are summarized in the upper right 
plot of Figure ^] which represents contours of constant 
star formation rates at z = and z = 3 in the r/ w - 
plane. The observed central value is shown as the bold 
lines. The two other contours correspond to a SFR value 
a factor of 2 higher and lower than the central value. As 
we see, the main constraint given by the observed SFR is 
that t„ must be lower than 5 Gyr. If not, the SFR will 
be underestimate at high redshift by a factor more than 
2. The contours also indicate that winds are required to 
reproduce the low SFR observed at low redshift. More 
qualitatively, winds are also required to prevent from 'the 
overcooling problem' and to remove baryons from the cold 
gas phase. 

5.3. Cosmic Stellar Density 

A complementary method to investigate the history of the 
baryon assembly in the universe is to observe the evo- 
lution of the global stellar mass density £l*(z). This is 
an independent way of constraining the model, since the 
SFR is related to young stars, while cosmic stellar obser- 
vations focus on old, red and low mass stars. The middle 
left plot of Figure fT 5| shows the obse r vation al data points 
of fi» compiled bv IDickinson et all ll2003l) from various 
near-infrared and optical observations et al. 2001; 

Brin chmann & Ellia l2000t ICohenl |2002; Di ckinson et al.l 
20031) . 

In order to perform an accurate comparison, we need 
to take into account the death of short-lived stars during 



the history of the universe. Since most of the stars are 
formed when the universe is a few Gyr young, this effect 
is likely to be important. We use h ere the stellar Initia l 
Mass Function (IMF) proposed by iKroupa et alJ {l993) 
in order to estimate the amount of stars still alive at a 
given redshift. The corresponding predicted stellar den- 
sity is plotted in the middle left plot of Figure ^] Our 
best fit model for this observations is now — 10 Gyr 
and 7/ w = 0.5, in complete disagreement with our best-fit 
model for the Cosmic SFR. This rather surprising result is 
due entirely to the high redshift observations of O* . Low 
redshift observations, on the other hand, are compatible 
with the SFR constraints. 

This p uzzle was identified as a possible 'missing galaxy 
problem' l)Nagamine et alJl2004h . This indicates that, if 
the IMF remains universal, the star formation rate should 
fall off at high redshift more strongly than in Figure El ln 
the model, this can be obtained for = 10 Gyr, but then, 
as can be seen in Figure IT7I the SFR is too high at low 
redshift. The possible inconsistency between the observed 
SFR and the observed ste l lar de nsity is discussed in great 
details in lDickinson et all l)2003|) . The mass- luminosity re- 
lation might be poorly estimated at high redshift because 
the galaxy luminosity function is not well constrained. 
Another solution is to invoke an evolving IMF, with more 
high mass stars at high redshift. Finally, the last solution 
is to invoke a new physical process, not included in the 
model, which inhibits star formation at high redshift. 

5.4. Extragalactic Background Light 

Another strong observational constraint is the 
Extragalactic Background Light (EBL), integrated 
from the UV to the IR. The EBL is an estimate of the 
total amount of energy emitted by stars and AGN in 
the history of the universe. Consequently, the cumulative 
luminosity of all stars created b y the model cannot exceed 
the value of this background. IMadau fc Pozzettil l|200f]h 
compute the value of the observed EBL, integrating 
from 0.2 to 2000 [irti, and found Irb t , = 55 nW/m 2 /sr 



Following the method presented in Madau fc Pozzettil 
(2000?), we compute the EBL corresponding to our 
canonical model 



Iebl — 



47T 



Pbo\(t) 
l + Z 



dz, 



(41) 



with the bolometric emissivity at epoch t 

pbol(t)= [ L(T)p*{t~T)dT. (42) 

Jo 

In this expression, L(r) is the bolometric luminosity of a 
single stellar cluster of unit mass, as a function of the 
cluster age r. W e use for L(t) the analyt ical approx- 
imation given by IMadau k. Pozzettil l|2000|) for a solar 
metallicity stellar population. We find for our canonical 
model /ebl — 100 nW/m 2 /sr, larger than the expected 
value, but within ob servational uncer tainties (for exam- 
ple optical data from IBernsteini l|l999h and IR data from 
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Fig. 18. Upper left plot: Cosmic Sta r Formation Rate a s a function of reds h ift. Data po i nts, uniformally corrected from 
various observational biases, come f romlHughes et al.1 Il99sj):ISteidel et~all Jl999l) : lFlores et alJ (ll999h : lGlazebrook et alJ (119991) : 
lYan et all lll999l) : lMassarotti et alJ (l2001^ : lGi£wahscoet^dT 1120041) . The solid line is our canonical model (i* = 3 Gyr, r? w = 1.5). 
The dashed line is our best fit model for the SFR (f* = 1.5 Gyr, 7/ w = 1). Upper right plot: contour of constant star formation 
rate (solid contours: z — 0, dashed contours: z = 3) as a function of t, and r/ v . Bold contours are the observed value, while 
others correspond to a factor of 2 above and below the observed centr al value. The cross is th e position of our canonical model. 
Middle left plot: Cosmic Stellar Density as a function of redshift, from lDickinson et al.l ll2003f) . The solid line is the prediction of 
our canonical model, while the dashed line is our best fit model (£» = 10 Gyr, rfc = 0.5). Middle right plot: contour of constant 
stellar densit y (solid contours: z = , dashed contours: z = 3). Lower left plot: cosmic gas density in damped Lyman alpha 
systems, from lSomerville et alJ i200ll) . The solid and dashed lines corresponds to the mass fraction in cold gas, as predicted by 
our canonical model. Lower right plot: contour of constant cold gas density (solid contours: z = 0, dashed contours: z = 3). 
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Fig. 19. Top: average star form ation history as a function of lookback time for different final stellar masses. The observed 
values from lHeavens et alJ i2004) are plotted on the top left panel. As in the article, the curves are offset vertically successively 
by 0.5 in log except for the most massive galaxies which are offset by an additional 1.0. The top right panel shows the average 
star formation history of our canonical model (t* = 3 Gyr and ry w = 1.5). For the most massive galaxies, a second model with 
metal-rich cooling and 'superwinds' is also plotted, in better agreement with the data. Bottom: mass fraction in the hot X-ray 
emitting gas in several observed groups and clusters (ISanderson et alJl2003h . The lower left solid line is the prediction of our 
canonical model, while the lower right solid line is the prediction of our 'superwind' model. 



ICharv & Elba3 (|200l[) and iFixsen et all {1998) lead to 
Iebl — 80 nW/m 2 /sr). Our best fit model for the Cosmic 
SFR (t* ~ 1.5 Gyr and t] w ~ 1) gives Jebl — 150 
nW/m 2 /sr, and can be therefore considered as being ruled 
out by the observational constraint. 

5.5. Cosmic Baryon Budget 

T he observed Baryon Bu dget is discussed in great details 
in iFukueita et all ([j.998). They infer from various obser- 
vations at z = baryon mass fractions of fiback — 0.002, 
ft» = 0.0035, f2 co id = 0.00063 and fl hot = 0.017. We note 



immediat ely that, if one consi ders the most recent WMAP 
estimate l|Spergel et al.ll2003|) of O b ~ 0.04, roughly 50% 
of the baryons are missing. The most striking disagree- 
ment between these estimates and our model is the very 
low baryon fraction in the background (or Lyman alpha 
forest) at z = 0. As it is now admitted, Lyman alpha 
observations at low redshift are very uncertain, so we 
consider here that most of the missing baryons are in 
fact in the diffuse b ackground. Recent observations from 
IPenton et alJ l)2004() seems to confirm that at least 30% of 
the baryons lie in the Lyman alpha forest. 
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Another conseque nce of the analysis performed by 
iFukugita et al.1 1^998) is that baryons in the condensed 
phase (cold gas + stars) are only 10% of the total amount 
of baryons in the universe (fi b — 0.04). This strongly sup- 
ports the presence of galactic winds, in order to overcome 
the 'overcooling problem'. Additional support is given by 
the fact that a large fraction (~ 40%) of baryons are today 
in the hot gas phase that is to say plasma in groups and 
clusters. Our canonical model (with winds) predicts frac- 
tions (relative to / b ) of / bac k ^ 50%, / ho t ^ 30%, f ro id ^ 
1% ar id /* ~ 20% (10% with the IMF of lKroupa et all 
l)l993|) L in rough agreement with present day observa- 
tions. Without winds, the same model is much more diffi- 
cult to reconcile with observations, since we obtain in this 
case / back ~ 40%, / hot ~ 20%, / cold ~ 4% and /, ~ 40%. 

Another strong observational constraint come from the 
evolution of the cold gas density, deduced from the obser- 
vations of Damped Lyman Alph a, Systems (DLAS ) indis- 
tant quasars spectra. Following ISomerville et al.1 ll200ll) 
we use obse rvations perfo r med b v lStorrie-Lombardi et all 
l)l996|) and IZwaan et alJ l|l997|) to estimate the evolu- 
tion of O co id as a function of redshift. These observa- 
tions of DLAS are lo wer limit estimates of O co id (see 
ISomerville et alJ l)200lh for a discussion). Dust absorption 
is likely to have a strong effect in biasing these results. 
We use the method proposed by IPei et all l|l999|) to cor- 
rect from these effects. 

Observational data points are shown in the lower left 
plot of FigureEl The curve reaches its peak value O co id — 
0.004 at a redshift z ~ 2, in good agreement with our 
canonical model, which, in this case, is also the best fit 
model. It is worth noticing that the observed f2 co id curve 
is proportional to the observed SFR curve. This can be 
considered as a nice consistency check in the observational 
data and provides support to our simple star formation 
model. 

The lower right plot of Figure El shows contours of 
constant il co id at z = and z = 3 in our model parameter 
space. Contours appear as straight lines in the f*-ry w plane. 
This is consistent with the fact that the cold gas fraction 
depends mainly on the gas depletion time scale which can 
be estimated to be i* /(??«, + 1) (see Section l3~7)) . The low 
observed value of £l co id favors a small depletion time scale 
with 1 < t*/(i} w + 1) < 3 Gyr, a rather tight constraint. 

5.6. Halo Star Formation History 

We now analyze recent observations of individual galaxies 
star formation history. Using Extended Press & Schechter 
theory, we can apply our model to predict the baryon his- 
tory within individual halos, in an a verage sense. Us i ng the 
Sloan Digital Sky Survey (SDSS), iHeavens et all <j2004^ 
infer from 10 5 nearby galaxy optical spectra the age distri- 
bution of their stellar population. For each galaxy, they de- 
duce its individual star formation history (using a Salpeter 
initial mass function). Finally, they compute the average 
SFR of all galaxies in a given stellar mass range. This 



quantity can be directly compared to the prediction of 
the model, if one converts the halo Virial mass into a halo 
star mass, using our model 'star-to-mass ratio'. 

Interestingly, the observed average star formation his- 
tories (upper left plot of Figure IH^I seem to indicate that 
large galaxies (M* > 10 12 Af Q ) form stars earlier than 
small ones (M* < 1 10 Mp>). This beha vior was called 
'anti-hierarchical' bv lHeavens et all l)2004|) and was iden- 
tified as a potential problem for the hierarchical scenario 
of structure formation. The upper right plot of Figure El 
shows our canonical model predictions (t* — 3 Gyr and 
rjvj = 1.5). The exact quantity we plot is 

P* = U {M , z )M n(M )AM, (43) 

where n(Mo) is the Press & Schechter halo mass function. 
We see that the same 'anti-hierarchical' behavior (down- 
sizing effect) is reproduced by our model, within the hier- 
archical collapse framework . There are two explanations. 
First, low mass halos have a total mass Mq very close to 
the Minimal Mass Af m i n , so that the mass fraction in 'star 
forming progenitors' is smaller than for high mass galaxies. 
This can be seen in Figure [5J where the diffuse gas mass 
accretion rate vanishes as M — > M min . Second, for large 
mass halos, the cooling efficiency drops at low redshifts, 
as more and more progenitors reach T max , the maximum 
cooling temperature. At this point, no more fresh gas is 
accreted onto the cold disc and star formation is slowing 
down. 

The observed individual star formation histories are 
therefore successfully reproduced by the model on a qual- 
itative level. If one looks carefully into Figure El one sees 
that for large galaxies, the predicted star formation history 
disagree with the observed one. The sudden drop in the 
star formation rate (interpreted as the end of gas cooling 
into cold discs) happens earlier in the model (look-back 
time between 7 and 8 Gyr) than in the data (look-back 
time around 2 Gyr). Moreover, after gas cooling ends, the 
decrease in the star formation rate is much faster in the 
data than in the model. 

The first feature can be accounted for if one takes 
into account metal enrichment into the cooling gas. For 
a metallicity of one tenth or one third solar, cooling can 
be significantly higher than for a primordial H and He 
plasma. We have tried a new model with an increased 
Maximum Cooling temperature T max = 2 x 10 6 K. In 
this case, the position of the knee in the star formation 
rate curve agrees perfectly with observational data. On 
the other hand, the strong decrease after the knee is dif- 
ficult to explain. Note however that the measure of the 
SFH from the optical spectrum is a very complex opera- 
tion and the uncertainties are important. For example, a 
large part of the star formation is enshrouded and may be 
missed. Moreover, the IMF is probably not a Salpeter one. 
Consequently, the decrease might be less accentuated than 
presented. Nevertheless, if the decrease is confirmed, one 
solution is to invoke very strong winds that remove most 
of the cold gas accumulated before the end of disc accre- 
tion. To illustrate this, we further modify the canonical 
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model, introducing what we call here 'superwinds', with 
parameters T w ~ 10 7 K and rj w = 15. This last model 
(see Fig. I19|) is finally able to reproduce the observed star 
formation history in large galaxies. 

This 'superwinds' scenario is completely ruled out by 
the observed global baryon history (see previous section). 
We have therefore to consider 2 coexisting wind models 
of very different nature: 'galactic winds', driven by super- 
novae bubbles, for normal and dwarf galaxies, and 'su- 
perwinds', possibly driven by a massive central black hole 
or AGN for massive galaxies. Building a self-consistent 
model along those lines is beyond the scope of this paper, 
but several attempts of 's uperwinds' models c an already 
be found in the literature IjSpringel et al. 2005). 

5. 7. Hot Gas Fraction in X-ray Clusters 

The observed amo unt of hot gas i n grou ps and clusters 
(fihot = 0.017 in iFukugita etafl l|l998|U put a strong 
constraint on the hierarchical scenario, mostly by requir- 
ing galactic winds to solve the 'overcooling problem'. 
ISanderson et al.l (j20Q3^ have analyzed various X-ray ob- 
servations of groups and clusters and computed the frac- 
tion of hot gas (/hot) as a function of the observed X-ray, 
emission weighted, temperature. These data points are 
plotted with their corresponding error bars in Figure El 

First, we learn from these observations that most 
baryons in clusters have to be in the hot phase, unless 
we are ready to face a serious crisis with the WMAP con- 
straint Sib = 0.04. Second, there is a clear correlation be- 
tween the fraction of hot gas in each halo and the X-ray 
temperature. Moreover, this correlation /hot oc Tx is often 
used as a possible explanation for the observe d Lx — Tx 
relation in clusters ((Neumann fc Arnaudll20oH) . 

We have plotted on the same figure our canonical 
model prediction for the hot gas fraction as a function 
of the Virial temperature. In small galaxies, where both 
galactic winds and gas cooling are important, the hot 
gas fraction have its minimum around /hot — 3%. For 
large halos, cooling is less and less efficient, and in the 
same time, winds can not escape the halo potential well. 
This double mechanism (cooling + winds) has the im- 
portant consequence of refilling with hot gas the parent 
halo. This qualitative picture is interesting, but when one 
compares our canonical model predictions with the X-ray 
data, the result is quite disappointing. The refilling mecha- 
nism we have just explained occurs at too low temperature 
(T ~ 0.1 keV), while data suggests a significantly higher 
transition temperature. Note however that the observa- 
tions concern the center of the clusters and the extrapola- 
tion t o larger radii (using [3 model) is not safe IjNeumannl 

One solution might come again from our 'superwinds' 
scenario. We modify our canonical model by first increas- 
ing the Maximum Cooling temperature up to T max = 
2 x 10 6 K, as it should be for the case of a realistic, 
metal-rich, plasma. We then modify our wind parameters 



T 




M (h" 2 M ) 

Fig. 20. Predicted and observed stellar and HI mass function. 
The thick dashed line is the stellar mass function predicted by 
our standard model, whereas the thin dashed line is for our 
superwind model. Shown as dia monds is the Sch echter fit to the 
observed stellar mass function JCole et all200lL Similarly, the 
thick dot-dashed line is the cold gas mass function predicted 
by our standard model whereas the thin dot-dashed line is for 
our s uperwind model. Th e Schechter fit to the HI mass function 
from Zwa an et al] (J2005) is shown as triangles. For comparison, 
the thick continuous line is the Press-Schechter mass function, 
assuming a universal baryon fraction everywhere. 

to T w ~ 10 7 K and r/ w = 15, as the 'superwinds' model 
we use in the last section. We plot the hot gas fraction 
we obtain for this new model in Figure IHH The transition 
from 'gas poor' to 'gas rich' regime occurs now at a much 
more realistic temperature (T ~ 1 keV). This transition is 
however still much sharper in our model than it is in the 
data. As suggested in the last section, we could therefore 
improve our wind model by explicitly introduce 2 different 
fee dback scenarios: ' supernovae driven' and 'AGN driven'. 

iKavet al.lll2004h have performed hydrodynamical sim- 
ulations of galaxy clusters with a kind of feedback which 
is very close to our superwind scenario. Indeed, they heat 
the dense and cold gas of their clusters at a temperature of 
17 keV. This corresponds in our notations to T w ~ 2 x 10 7 
K (compared to T w ~ 10 7 K in our superwind model). 
Using this strong feedback, they reproduce the observed 
cluster Lx — Tx relation with the correct level of entropy 
in the ICM core. Such a strong feedback seems therefore 
essential to reproduce both global fraction of hot gas and 
entropy profile. 

5.8. Stellar and HI mass functions 

In this final part, we investigate the stellar and cold 
gas mass function at z=0. We reach here the limit of 
our model because, as mentionned before, our analyti- 
cal predictions are for baryon mass fractions as a func- 
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tion of halo masses, whereas observed baryon mass frac- 
tions are usually given as a function of individual galaxy 
masses. The conversion between galaxy and halo masses 
can be performed usi ng the 'halo occupation number' 
l|Kravtsov et alJl2004aj) . which gives the average number 
of galaxy satellites within i?200 a s a function of the par- 
ent halo mass. We intend to apply this correction to our 
model predictions in a future paper. Nevertheless, we com- 
pare our predicted stella r and cold gas ha lo mass function 
to the observe d stellar llCole et al.ll2fl0l|) and HI galaxy 
mass function l|Zwaan et alJl2005|) 7 

As presented in figure I2U1 the global normalisation of 
both stellar and cold gas mass function is in good agree- 
ment with the observations. Indeed the integral of this 
mass function multiplied by the mass is nothing else that 
the z=0 cosmic baryon budget for this 2 phases. As for 
the observations, the profile shape shows a fall-off for high 
mass and a flatening for lower masses. This is easily ex- 
plained by our predicted stellar and cold gas mass fraction. 
For the high mass tail, the fall off is due to the fast decline 
of Press-Schechter mass function, where as for low mass 
the flatening is due to the minimal halo mass M m j n above 
which the fraction of baryons become weaker and weaker. 

If we look in more details, we see an encouraging agree- 
ment for intermediate mass. The decline of the high mass 
tail is however steeper in the observations, both for the 
HI mass function r(Mhi) and the stellar one n(M*). We 
have identified 2 reasons for that. First, as highlighted 
before, the halo occupation number should be taken into 
account, because it would decompose each large halos in 
a collection of lower mass halos. As a consequence, the 
high-mass decline may be more steep. Second, superwinds 
may also lo wer the cold and ste llar mass fraction in high 
mass halos ijSpringel et al.ll2005|) . To illustrate this point, 
we have plotted the effects of our superwind scenario on 
the mass function. Note that a more complex model with 2 
different winds would result in more realistic predictions. 
For low mass galaxies, the model tends to underestimate 
the cold gas mass function. Here again the halo occupa- 
tion number would explain part of this discrepancy by 
increasing the number of low mass halos. In our model 
the transition between the star forming halos and the dif- 
fuse background occurs abruptly at M m i n . As seen in our 
simulations, this transition is in fact much smoother: this 
is likely to improve the agreement between the predicted 
and the observed mass function. 

6. Conclusion 

We have studied the baryon budget evolution in the frame- 
work of the hierarchical scenario of structure formation, 
using 4 different phases: diffuse background (Lya Forest), 
hot gas (plasma in the halo of galaxies, groups and clus- 
ter), cold discs and stars. We have paid particular atten- 
tion to the star formation rate, as it is a key observational 
constraint for our current cosmological model. 

We have analyzed the baryon history for the universe 
as a whole, but also on a halo-by-halo basis. For that pur- 



pose, we have developed a fully self-consistent (though 
simple) analytical model. These last two point are the 
most original aspects of our work. Our analytical model 
has proven to be an efficient tool to quickly compute ac- 
curate predictions for the baryon budget history. It is cur- 
rently available as a set of IDL routines, and can be pro- 
vided by the authors upon request. 

In order to validate this model, we have performed nu- 
merical simulations of galaxy formation using the AMR 
code RAMSES. Our highest resolution run reach 512 3 
dark matter particles and half a billion AMR cells, which 
is among the largest galaxy formation simulations per- 
formed so far. We have also a nalyzed the simulation re- 
sults of ISm-ingel fc Hernauistl l|2003bh based on the SPH 
code GADGET. We found in all cases a good agreement 
between simulations and our model. This cross-validation 
has allowed us to use our model to analyze observational 
data. 

We have explored our physical parameter space i* — 
7/ w (star formation time scale and wind efficiency) and 
compared the model results to the cosmic observations 
of the comoving star formation rate, the evolution of the 
comoving density of stars and cold gas, and the intensity of 
the integrated extragalactic background. The conclusion 
is that the parameters i* = M co id/ = 3 Gyr and rj w = 
M^nd/ M* — 1.5 are favored. It means that winds with 
ejection rates around 1 or 2 times the star formation rate 
are required to prevent the overcooling problem. 

Comparisons with individual halo properties, such as 
the age distribution of stars in galaxies and the hot gas 
fraction in clusters seems to indicate that high velocity 
and high intensity outflows ('superwinds') are required in 
massive galaxies. The origi n of these violent out flows could 
come from a central AGN l|Springel et al.l2005|) . The mod- 
eling of such winds and their exact role in the metal en- 
richment of the IGM are currently under investigation. 
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Appendix A: Numerical methods 

A.l. Hybrid scheme for high-Mach flow 

Our hydrodynamical solver is built using a multidimen- 
sional unsplit Godunov scheme, with second order accu- 
racy both in space and in time. In order to capture cor- 
rectly shock waves, the total fluid energy is used as inde- 
pendent variable. As soon as velocity gradients (velocity 
undivided differences in each direction) in the flow have a 
comparable magnitude to the fluid sound speed, the fluid 
temperature can be computed accurately. In the oppo- 
site case, which happens routinely in cosmological appli- 
cations, the fluid temperature is dominated by truncation 
errors in the velocity field. Since the velocity field is mainly 
determined by the gravitational acceleration, this arises 
when the mesh resolution Ax is greater than the Jeans 
length. In order to deal with this problem, particularly 
acute in our case, since we need an accurate temperature 
determination for the cooling and heating functions, we 
are left with two possibilities 

— we refine our mesh sufficiently in order to resolve the 
Jeans length. This solution is act ually adopted in most 
star formation numerical studies ( Truclovc et al 1ll997t 
iBoss et al~ll2000() . Except for our smallest boxes, this 
strategy is not applicable here. 

— we switch off the total energy conservative update in 
regions where truncation errors are too large, and use 
instead the internal fluid energy as independent vari- 
able. The problem is now to design the correct switch- 
ing strategy. 

Note that this problem, rather specific to cosmol- 
ogy, has been identifi ed by several auth ors, both in the 
Godunov fr amework llBrvan et al.lll995|) and in the SPH 
framework ( Shapiro et alJll99fit) . The solution adopted by 



c s — 7(7 ~ l) e cons to an estimate of the local velocity 
truncation errors Axd x u 



these various authors was to carefully detect shock waves 
formation, in order to turn off spurious heating (or cool- 
ing) before shock heating actually occurs. In RAMSES, 
we have adopted the following approach. At the end of 
the hydro step, we update the fluid internal energy e = pe 
by solving in parallel two different sets of equations. 

d t E + d x u(E + p ) = and e = E-^pu 2 , (A.l) 



<9 t e + d x ue = —Pd x it = — (7 — l)ed x u. (A. 2) 

The first one is the standard conservative update, followed 
by the subtraction of the fluid kinetic energy, while the sec- 
ond one is a direct, non conservative update of the internal 
energy. This hybrid approach is a classical trick in fluid 
dynamics, to overcome certain problems a ssociated to the 
conservative formulation in equation IA.1I ( Torolll99"?0 . In 
order to solve for eauation lA.2l we need to compute (and 
store) the velocity divergence d x u for each fluid element. 
We then decide which update between e pr i m and e cons is 
the valid one by comparing the conservative sound speed 



Cfinal — ^cons 
^final — Cp r i m 



if 



> a{Axd x u) 2 



(A.3) 



otherwise. 



a is a free dimensionless parameter of order unity. We take 
as phase value a — 0.5 in all the simulations presented 
he re. Note that our ap proach differs from the one proposed 



i lBrvan et al.l 



in several aspects, the most impor- 
tant one to our opinion is that criterion IA.3I is Galilean 
invariant as for the Euler equations. 

The consequence of our hybrid energy update is that 
shock heating is turned off completely in small mass ha- 
los, where truncation errors dominate the fluid internal 
energy, due to a lack of resolution. This introduces an 
effective low mass cut off in the simulated halos distribu- 
tion, below which the temperature evolution is governed 
only by adiabatic compressional heating and UV heating: 
small mass halos are therefore considered as sitting in the 
smooth, Lyman alpha background. Careful numerical ex- 
periments with RAMSES, under the refinement parame- 
ters used in this paper, indicate that this minimal mass 
sits around 100-200 dark matter particles. 

A. 2. Refinement strategy 

We adopt here a refinement strategy based on the clas- 
sical 'quasi-lagrangian' approach. Our goal is to describe 
as accurately as possible the collapse of individual Fourier 
modes. We therefore refine cells when the mass contained 
in that cell exceeds a certain threshold. When only one 
fluid is present (dark matter say), this threshold trans- 
lates into a critical particle number, above which new 
refinements are triggered. This approach has been used 
by several author s performing AMR N body simulations 
in the literatu re l|Kravtsov et al.lll997r. libel et al.ll200d: 
iTevssierlEool. 

In the present case, we are dealing with 3 different 
fluids: dark matter, gas and stars. Moreover, 2 fluids are 
treated using particles (dark matter and stars), and one 
is described using grid cells. When cooling is included in 
the gas energy equation, baryons sink very efficiently in 
the center of their host halos, and ultimately dominates 
the total mass distribution. We have adopted the following 
approach in all the simulations performed here: we define 
the total number of 'particles' as 

AT tot = Ax 3 (p d m/m dm + p g /m g + p*/m*), (A.4) 

where TOdm is the mass of each dark matter particle, m g — 
/bTOdm is the initial gas mass of each cell in our base grid 
(,/b = f^b/fim), and to* is the mass of each star particle. 
We set to* = TO g , to prevent spurious refinements or de- 
rcfincment in star forming regions. In this way, the grid, 
as well as the source term for the Poisson equation, are 
not aware of the underlying star formation process. 

AMR cells are refined, at each level of refinement, 
if the effective number exceed AT to t = 40. This thresh- 
old is rather high compared to previous N body AMR 
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simulations where TVt.nt ^ 5 — 10 l|Kravtsov et al.lll997t 
lAbel et alJbood: iTevssiedEx)^ . In our current gas and 
dark matter runs, we want to minimize the appearance of 
gas cells devoid of dark matter particles. This is known to 
lead to spurious oscillations in the gas density l|Tevssierl 
l2f)02j) . 

We set the maximum level of refinement to ^ max = 5. 
This number is rather low compared to p revious N body 
AMR simulations, where l ma , y — 7 — 8 l|Kravtsov et al.l 
ll997tlAbel et a.l]l200(i fTevssierll2002|) . Since we are deal- 
ing with strong cooling in high density regions, there is 
nothing to prevent the gas from collapsing into dense cold 
knots. The AMR grid will be automatically refined, lead- 
ing to very small cells. If our maximum refinement level is 
too high, dark matter particles sitting in gas dominated, 
high density regions will suffer from two body relaxation 
effects. Setting £ max = 5 gives us a good compromise be- 
tween spatial resolution and cautiousness. 

As one can see, our refinement strategy is rather cau- 
tious, and based on a very conservative choice of parame- 
ters, compared to previous, more aggressive strategies pro- 
posed in the AMR literature. For our largest runs, start- 
ing with a base grid of 512 3 coarse cells, this allows us 
to reach a formal resolution of 16384 3 in the densest re- 
gions of the flows. This configuration was computed with 
256 processors and 350h wall-clock time on a massively 
parallel computer. 

Appendix B: Filtering Mass 

Thoul & Weinberg (97) and lGnedinl l(2000t were the first 
authors to carefully examine how baryons fall (or not) 
into dark matter halos potential wells in presence of a 
UV background. These authors realized that it is only 
above a critical mass, called the 'Filtering Mass' Mp, that 
baryons and dark matter assemble in comparable amount 
inside virialized halos. In halos of mass lower than Mp, 
the baryon fraction was found to be (on average) lower 
than the universal value. This means that below this mass, 
baryons can be considered as a quasi-homogeneous, diffuse 
co mponent . 

ICnedinl l(2000h noticed that the instantaneous Jeans 
mass Mj was not the correct way of computing Mp. The 
correct estimate of Mp at a given epoch is based on the 
average of the Jeans mass over the past history of the uni- 
verse. We therefore need a model for the thermal history of 
the background, in other words the average temperature 
f. 

At very early epochs, due to residual electrons, the 
gas temperature is tightly coupled to the cosmic black 
body. After z ~ 200, the gas temperature finally decou- 
ples and evolves as (1+z) 2 . At a given redshift (called here 
the reionization redshift) z r , we assume that baryons are 
promptly heated to a fixed temperature T T ~ 10 4 K, and 
remains roughly isothermal afterwards. Note that a more 
realistic computation of the ionization fraction would lead 
to a more complex behavior for the temperature after 
reionization. Our model for the background temperature 



is therefore given by the following equations 

( 2.73(1 + z) z > 200, 
f= I 0.0136(1 + z) 2 200 >z>z r , (B.l) 
I T, z < z T . 

The two parameters z r and T r must be carefully chosen in 
order to recover the proper therm al history of the universe . 
Our simulations are based on the lHaardt k Madaul l|l996f) 
model for the UV background. This rather complex model 
can be roughly reproduced with T r ~ 6 x 10 3 K and z r ~ 6. 
For a different model with a harder UV spectrum, the 
corresponding reionization temperature has to be set to a 
higher value. Recent observations by the WMAP satellite 
suggest a larger value of the reionization redshift z r ~ 20. 
This latter value will be used as a reference parameter for 
comparison to observational data. 

As soon as the background temperature is known, we 
can compute the Jeans mass 

where the Jeans scale is given by 

and p is the average mass density in the universe. Finally, 
the Filtering Mass is simpl y given by the following integral 
over the expansion factor l)Gnedinl l2000') 

M} = - f da'Mj (a')(l ~ (B.4) 
a Jq V a 

Thus, dark matter halos with masses lower than Mp (or 
equivalently with Virial temperature lower than Tp) are 
tracers of the diffuse baryon background, while dark mat- 
ter halos with masses greater than Mp contain an hot 
baryon component in hydrostatic equilibrium, whose mass 
fraction is close to the universal one /b = fib/^m — 0.13. 

Appendix C: Cooling model 

The Filtering Mass is the first important mass scale. It is 
however not the only one needed to properly define star 
forming halos. We need to determine a second mass scale 
above which hot gas in hydrostatic equilibrium is able to 
cool and fall in the center of the host dark matter halo. 

If cooling is ignored, this hydrostatic plasma can 
indeed be considered a s another diffuse component. 
Numerical simulations llEke et alJ Il998l). theoretical 
arguments llSuto et alJ ll99St iKomatsu k Seliakl 1200 it 
lAscasibar et al.ll2003l) and observations of X-ray e mission 
in large X-ray clusters ijNeumann k Arnaudll999|) all sug- 
gest that the maximum overdensity during adiabatic col- 
lapse lies around 10 , which is a rather low value compared 
to typical overdensities of galactic discs. The diffuse back- 
ground is therefore defined as all dark matter halos with 
masses lower than the Filtering Mass, plus all dark matter 
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Fig. C.l. Ratio of free-fall time over cooling time for a typical 
core gas density in the t emperature and redshift plane, for the 
lHaardt fc Madaul (1996) rcionization model (z r ~ 6). The up- 
per solid line shows the redshift evolution of the maximum cool- 
ing temperature T max , while the lower solid line corresponds to 
the minimum cooling temperature T coo i. The region enclosed 
by these 2 solid lines has t coo i < tg. Also shown as the dashed 
line is the evolution of the Filtering Mass, in units of Virial 
temperature. 

halos whose hot plasma is not able to cool and condense 
as a disc. 

In this paper, we only consider atomic cooling pro- 
cesses. We ignore molecular cooling processes. There are 
two reasons for tha t. First, we do n ot c onsider Pop III 
star formation, as in lAbel et alJ <)2002f) and lSokasian et all 
(200^), since it is currently a matter of debate to deter- 
mine if this first generation of stars has a negative or a pos- 
itive feedback on molecular hydrogen formation. Second, 
in this paper, we perform numerical simulations that do 
not have enough resolution to describe properly molecular 
cooling and fragm entation inside HI discs , as opposed to 
the recent work of iKravtsov et alJ l|2004a|) . 

Semi-analytical models of galaxy formation have de- 
veloped various methods to compute the atomic cooling 
efficiency as a function of halo mass and merging history. 
The cooling model we use here is based on a major sim- 
plification of such models, in order to deal with analytical 
formulae. We take the typical core density in the halo 
gas distribution to be roughly constant with S COIC ~ 10 5 . 
Using cooling; and heati ng functions corresponding to the 
lHaardt & Madaul l|l99fi() UV background, we compute the 
ratio of the net cooling time over the local free-fall time 
t C ooi /t ft i assuming isothermality at the Virial temperature 
of the halo. 

We plot in Figure IC~T1 this ratio as a function of red- 
shift and halo temperature. We also plot the evolution of 
the Filtering Mass (in units of Virial temperature) as the 



dashed line. Darker contours are for halos with very fast 
cooling rates, while lighter ones are for halos with cooling 
time slower than the dynamical time. 

The transition between fast and slow cooling occurs in 
this plot at the curve defined by t coo \ = tg. This curve 
can be approximated roughly by 2 straight lines, with a 
maximum cooling temperature given by 

T max ~6x 10 5 (l + z) 3 K, (C.l) 

and a minimum cooling temperature given by 

T cool ~ 6 x 10 3 K. (C.2) 

We finally compute the accretion rate of hot gas into the 
central, cold, centrifugally supported disc, for these 2 dif- 
ferent regimes 

— Fast cooling: cooling is assumed to fast enough so that 
mass accretion onto the disc is purely limited by the 
orbital decay of infalling gas clumps. For purely radial 
orbits, the time scale associated to this orbital decay 
is equal to t OYD = i?2oo/^200- For more complex orbits, 
the orbital path is larger, and so is the associated time 
scale. We introduce a new free parameter, the 'mean 
orbital length', and compute the associated orbital de- 
cay time scale as t OIO = -Rorb/V^oo- The cooling rate is 
finally given by w coo i ^ l/*orb- 

— Slow cooling: cooling is the limiting process that con- 
trols mass accretion onto the disc. We simplify even 
more by discarding completely the cooling flow de- 
scription, and assume uj coo \ ~ 0. 

In this analytical model, we have introduced the free pa- 
rameter i? rb that allows us to vary the orbital time scale 
relative to purely radial accretion. In the latter case, the 
time scale reduces to i?20o/^200; which is independent of 
halo mass. In the following, we assume that i? or b is propor- 
tional to -R200, so that i rb can also be considered as being 
independent of halo mass. This unknown (but constant) 
ratio i? O rb/-R200 will be determined using our numerical 
simulations as calibration tools. Note that this simple 
model can be improved in many ways. Semi-analytical 
work are good examples of a more accurate modeling. 
However, our simple model offers the possibility to per- 
form analytical calculations and turns out to be reason- 
ably accurate. 
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